J'ai complété le document sur la mise en orbite d'un satellite et ajouté en particuli...
[pst-eqdf.git] / gravitation / gravitation_01.tex
index 93aaa0b..e364ae1 100644 (file)
-\documentclass{article}
-\usepackage[a4paper,margin=2cm]{geometry}
-\usepackage[latin1]{inputenc}
-\usepackage{pstricks-add,pst-eqdf,pst-func}
-\usepackage{array,amsmath}
-\usepackage{listings}
-\newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.125,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.5 0}}}
-%%%%%%%%%%%%%%%%%%
-\title{Interaction gravitationnelle avec PSTricks}
-\date{19 juin 2\,012}
-\begin{document}
-\maketitle
-\def\eqsatellite{%
-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)}%
-\section{Mise en orbite d'un satellite}
-\begin{figure}[htbp]
-  \begin{center}
-\psset{unit=2}
-\begin{pspicture}(-6,-2)(2,5)
-\pstVerb{
-    /GM 1 def
-    /theta0 -35 def
-    /r0 1 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 1.3 def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def
-    /Lc r0 v0 mul def % moment cinetique
-    /par Lc dup mul GM div def % paramètre de l'ellipse
-% excentricité
-    /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
-%%%%%%%%%%%%%%
-    /a_2 par 1 exc dup mul sub div def % demi-grand axe
-    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
-    /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def
-% vitesse à l'apogée
-    /vA GM par div sqrt 1 exc sub mul def
-% vitesse au périgée
-    /vP GM par div sqrt 1 exc add mul def
-% coordonnées de vA
-    /vAx vA theta0 90 add cos mul neg def
-    /vAy vA theta0 90 add sin mul neg def
-% coordonnées de vP
-    /vPx vP theta0 90 add cos mul def
-    /vPy vP theta0 90 add sin mul def
-% Apogée
-    /xA par 1 exc sub div theta0 cos mul neg def
-    /yA par 1 exc sub div theta0 sin mul neg def
-% Périgée
-   /xP par 1 exc add div theta0 cos mul def
-   /yP par 1 exc add div theta0 sin mul def
-% Centre
-   /xO xP xA add 2 div def
-   /yO yP yA add 2 div def
-}%
-\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
-\psdot[dotstyle=+](0,0)
-\uput[d](0,0){O}
-\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-2)(2,5)
-\rput(1,4){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{0000}s}}
-\parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
- radius t cos mul
- radius t sin mul}
-\pscircle*(!x0 y0){0.05}
-\pnode(!xP yP){P} % périgée
-\pnode(!xA yA){A} % Apogée
-\pnode(!xO yO){O} % centre
-\rput(!x0 y0){\psline[style=vecteurB]{->}(!v0x v0y)\uput[ur](!v0x v0y){$\overrightarrow{v_0}$}}
-\psline[linestyle=dashed](A)(P)
-\pscircle*(A){0.05}
-\uput[ul](A){$A$}
-\uput[dr](P){$P$}
-% position du satellite à un instant quelconque
-\pstVerb{/theta_i 170 def
- /radius par 1 exc theta_i theta0 sub cos mul add div def
- /xS radius theta_i cos mul def
- /yS radius theta_i sin mul def
- /ux theta_i cos 1 mul def
- /uy theta_i sin 1 mul def
- /xi xS ux sub def
- /yi yS uy sub def
- /xi2 xS ux 2 div sub def
- /yi2 yS uy 2 div sub def}%
-\pnode(!xi2 yi2){Mi2}
-\pnode(!xi yi){Mi}
-\pnode(!xS yS){S}
-\pscircle*(S){0.05}
-\psline[linestyle=dotted](S)(0,0)
-\psline[style=vecteurA]{->}(S)(Mi)
-\uput[l](S){S}
-\uput[u](Mi2){$\overrightarrow{F}$}
-\psarcn{->}(0,0){0.4}{0}{!theta0}
-\uput{0.5}[!theta0 2 div](0,0){$\theta_0$}
-\rput(A){\psline[style=vecteurB]{->}(!vAx vAy)}
-%\rput(P){\psline[style=vecteurB]{->}(!vPx vPy)}
-\psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,5)
-\uput[r](0,4.9){$y$}\uput[u](1.9,0){$x$}
-\psdot(O)\uput[r](O){$\Omega$}
-\rput{!90 theta0 add}(O){\psline[linestyle=dashed](!b_2 neg 0)(!b_2 0)}
-\end{pspicture}
-\caption{Mouvement d'un satellite}
-  \end{center}
-\end{figure}
-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 :
-\[
-\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}
-\]
-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 :
-\[
-r=\frac{p}{1+\mathrm{e}\cos(\theta-\theta_0)}
-\]
-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.
-\[
-\mathrm{e}=\sqrt{1+\frac{2\mathcal{E}L^2}{mK^2}}\qquad p=\frac{L^2}{mK}
-\]
-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 :
-\[
-L=mr_0v_0\qquad \mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2
-\]
-En remplaçant $L$ et $\mathcal{E}$, on obtient pour l'excentricité et le paramètre les expressions suivantes :
-\[
-\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)}
-\]
-\[
-p=\frac{v_0^2r_0^2}{\mathcal{G}M}
-\]
-On se limite au cas du mouvement elliptique, avec, en conséquence, la condition :
-\[
-\mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2 < 0
-\]
-Le demi-grand axe $a$, le demi-petit axe $b$ sont :
-\[
-a=\frac{p}{1-\mathrm{e}^2}\qquad b=\frac{p}{\sqrt{1-\mathrm{e}^2}}
-\]
-La période $T$ qui obéit à la troisième loi de Képler :
-\[
-T^2=\frac{4\pi^2a^3}{\mathcal{G}M}
-\]
-La vitesse, en un point de l'ellipse, se calcule par :
-\[
-v^2=\mathcal{G}M\Big( \frac{2}{r}-\frac{1}{a}\Big)
-\]
-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 :
-\[
-v_P=\sqrt{\frac{\mathcal{G}M}{p}}(1+\mathrm{e})\qquad v_A=\sqrt{\frac{\mathcal{G}M}{p}}(1-\mathrm{e})
-\]
-
-\section{L'étude avec PSTricks}
-\subsection{La trajectoire}
-\def\parametres{
-    /GM 1 def % 4e14 def % GxM
-    /x0 6.5e6 def % position initiale
-    /y0 0 def
-    /vx0 0 def    % vitesse initiale
-    /vy0 1e4 def
-}
-%  x0   y0   x'0   y'0
-% y[0] y[1] y[2]  y[3]
-\[
-\left\{
-\begin{array}{rcl}
-\ddot{x}&=&-\dfrac{GM}{r^3}x\\[1em]
-\ddot{y}&=&-\dfrac{GM}{r^3}y\\
-\end{array}
-\right.
-\label{eq1}
-\qquad r=\sqrt{x^2+y^2}
-\]
-
-On peut dessiner la trajectoire du satellite et de ses caractéristiques de deux façons :
-\begin{itemize}
-  \item par l'utilisation de \verb+\parametricplot+ ;
-  \item ou celle de \verb+\psplotDiffEqn+.
-\end{itemize}
-\verb+\parametricplot+ utilise l'expression exacte de l'équation de la trajectoire en coordonnées polaires :
-\begin{verbatim}
-   \parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{%
-     /radius par 1 exc t theta0 sub cos mul add div def
-      radius t cos mul
-      radius t sin mul}
-\end{verbatim}
-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}.
-\begin{verbatim}
-\pstVerb{
-    /GM 1 def
-    /theta0 -45 def
-    /r0 0.5 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 1.92 def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def
-    /Lc r0 v0 mul def % moment cinetique
-    /par Lc dup mul GM div def % paramètre de l'ellipse
-% excentricité
-    /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
-%demi-grand axe
-    /a_2 par 1 exc dup mul sub div def % demi-grand axe
-%demi-petit axe
-    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
-% période
-    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def
-}%
-\end{verbatim}
-\verb+\psplotDiffEqn+ utilise les équations différentielles du mouvement, en notation algébrique :
-\begin{verbatim}
-%  x0   y0   x'0   y'0
-% y[0] y[1] y[2]  y[3]
-\def\eqsatellite{%
-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)}
- \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}%
-\end{verbatim}
-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.
-%
-\begin{center}
-\begin{pspicture}(-12,-2)(4,10)
-\pstVerb{
-    /GM 1 def
-    /theta0 -45 def
-    /r0 0.5 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 1.92 def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def
-    /Lc r0 v0 mul def % moment cinetique
-    /par Lc dup mul GM div def % paramètre de l'ellipse
-% excentricité
-    /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
-%%%%%%%%%%%%%%
-    /a_2 par 1 exc dup mul sub div def % demi-grand axe
-    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
-    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def
-}%
-\psframe*[linecolor=white](-3,-0.2)(0,0.2)
-\pscircle[fillcolor=blue!50,fillstyle=solid](0,0){0.5}
-\psgrid[unit=2,subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-1)(2,5)
-\rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}
-  \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}%
-\parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
- radius t cos mul
- radius t sin mul}
-\psdot[unit=2,dotsize=0.12](!x0 y0)
-\rput(!x0 2 mul y0 2 mul){\psline[style=vecteurC]{->}(!v0x v0y)}
-\end{pspicture}
-\end{center}
-\subsection{La vitesse}
-\verb+\psplotDiffEqn+ permet de voir comment varie la vitesse sur l'ellipse :
-\begin{verbatim}
-%  x0   y0   x'0   y'0
-% y[0] y[1] y[2]  y[3]
- \psplotDiffEqn[xunit=0.2,yunit=5,%
-                plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,
-                linecolor=red,method=rk4,plotpoints=1000,
-                algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%
-\end{verbatim}
-\begin{center}
-\begin{pspicture}(0,-1)(10,7)
-\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
-\pstVerb{/GM 1 def
-    /theta0 -35 def
-    /r0 1 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 1.3 def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def}%
-\psplotDiffEqn[xunit=0.2,yunit=5,
-  plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,
-  ,linecolor=red,method=rk4,plotpoints=1000,algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%
- \multido{\i=1+1,\I=5+5}{9}{\uput[u](\i,0){\I}}
-\pnode(! 36.4 5 div 0){P}
-\psdot(P)\uput[d](P){Périgée}
-\psline[linestyle=dashed](P)(! 36.4 5 div 7)
-\pnode(! 36.4 10 div 0){A}
-\psdot(A)\uput[d](A){Apogée}
-\psline[linestyle=dashed](A)(! 36.4 10 div 7)
-\uput[l](0,6.5){$v$}
-\psline[arrowinset=0.1,arrowsize=0.2]{<->}(10,0)(0,0)(0,7)
-\uput[u](10,0){$t$(s)}
-\uput[ul](0,0){0}
-\end{pspicture}
-\end{center}
-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 :
-\[
-\overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}
-\]
-$\dot{\theta}$ et $\dot{r}$ s'obtiennent par les relations suivantes :
-\[
-\dot{\theta}=\frac{r_0v_0}{r^2}
-\]
-\[
-\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)
-\]
-La chaîne de calculs est la suivante : $\theta\Longrightarrow r\Longrightarrow \dot{\theta}\Longrightarrow \dot{r}\Longrightarrow \overrightarrow{v}$.
-
-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.
-
-On sauve successivement le tableau des positions et celui des vitesses.
-\begin{verbatim}
-\psequadiff[method=rk4,plotpoints=1000,
-            algebraic,
-            whichabs=0,whichord=1,
-            tabname=XiYi]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
-\end{verbatim}
-\begin{verbatim}
-\psequadiff[method=rk4,plotpoints=1000,
-            algebraic,
-            whichabs=2,whichord=3,
-            tabname=vxvy]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
-\end{verbatim}
-
-Pour ensuite dessiner la trajectoire et les vecteurs-vitesse.
-\begin{verbatim}
-%\listplot[unit=1]{vxvy aload pop}
-% on dessine la vitesse un point sur 100
-\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
-\multido{\i=0+100}{20}{%
-\pstVerb{/vX vxvy \i\space get def
-         /vY vxvy \i\space 1 add get def
-         /xi XiYi \i\space get def
-         /yi XiYi \i\space 1 add get def}%
-\rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}
-\end{verbatim}
-\begin{center}
-\begin{pspicture}(-10,-10)(6,7)
-\psset{unit=2}%
-\pstVerb{/GM 1 def
-    /theta0 30 def
-    /r0 2 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 0.85 def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def
-    /Lc r0 v0 mul def % moment cinetique
-    /par Lc dup mul GM div def % paramètre de l'ellipse
-% excentricité
-    /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
-%%%%%%%%%%%%%%
-    /a_2 par 1 exc dup mul sub div def % demi-grand axe
-    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
-    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def}%
-\rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}
-\psequadiff[method=rk4,
-            plotpoints=1000,
-            algebraic,
-            whichabs=0,
-            whichord=1,
-            tabname=XiYi
-%           ,saveData,filename=XiYi.dat
-]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
-\listplot{XiYi aload pop}
-\psequadiff[method=rk4,
-            plotpoints=1000,
-            algebraic,
-            whichabs=2,
-            whichord=3,
-            tabname=vxvy
-%           ,saveData,filename=vxvy.dat
-]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
-%\listplot[unit=1]{vxvy aload pop}
-% on dessine la vitesse un point sur 100
-\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
-\multido{\i=0+100}{20}{%
-\pstVerb{/vX vxvy \i\space get def
-         /vY vxvy \i\space 1 add get def
-         /xi XiYi \i\space get def
-         /yi XiYi \i\space 1 add get def}%
-\rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}
-\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(3,3)
-\end{pspicture}
-\end{center}
-
-
-
-\section{Mouvement circulaire}
- \begin{center}
-\psset{unit=1}
-\begin{pspicture}(-5,-5.5)(5,5.5)
-\pstVerb{
-    /GM 1 def
-    /theta0 60 def
-    /r0 5 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 GM r0 div sqrt def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def
-    /Lc r0 v0 mul def % moment cinetique
-    /par Lc dup mul GM div def % paramètre de l'ellipse
-% excentricité
-    /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
-%%%%%%%%%%%%%%
-    /a_2 par 1 exc dup mul sub div def % demi-grand axe
-    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
-    /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def
-% vitesse à l'apogée
-    /vA GM par div sqrt 1 exc sub mul def
-% vitesse au périgée
-    /vP GM par div sqrt 1 exc add mul def
-% coordonnées de vA
-    /vAx vA theta0 90 add cos mul neg def
-    /vAy vA theta0 90 add sin mul neg def
-% coordonnées de vP
-    /vPx vP theta0 90 add cos mul def
-    /vPy vP theta0 90 add sin mul def
-}%
-\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.75}
-\psdot[dotstyle=+](0,0)
-\uput[d](0,0){O}
-\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(5,5)
-\rput(0,-2){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{000000}s}}
-\parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
- radius t cos mul
- radius t sin mul}
-\pscircle*(!x0 y0){0.1}
-%\pnode(!par 1 exc add div theta0 cos mul par 1 exc add div theta0 sin mul){P} % périgée
-%\pnode(!par 1 exc sub div theta0 cos mul neg par 1 exc sub div theta0 sin mul neg){A} % Apogée
-\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}$}}
-\uput[ur](!x0 y0){$M_0$}
-\psline[linestyle=dashed](0,0)(!x0 y0)
-% position du satellite à un instant quelconque
-\pstVerb{/theta_i 170 def
- /radius par 1 exc theta_i theta0 sub cos mul add div def
- /xS radius theta_i cos mul def
- /yS radius theta_i sin mul def
- /ux theta_i cos 1 mul def
- /uy theta_i sin 1 mul def
- /xi xS ux sub def
- /yi yS uy sub def
- /xi2 xS ux 2 div sub def
- /yi2 yS uy 2 div sub def}%
-\pnode(!xi2 yi2){Mi2}
-\pnode(!xi yi){Mi}
-\pnode(!xS yS){S}
-\pscircle*(S){0.1}
-\psline[linestyle=dotted](S)(0,0)
-\psline[style=vecteurC]{->}(S)(Mi)
-\uput[l](S){S}
-\uput[u](Mi2){$\overrightarrow{F}$}
-\psarc{->}(0,0){1}{0}{!theta0}
-\uput{1.1}[!theta0 2 div](0,0){$\theta_0$}
-\psline[arrowinset=0.1,arrowsize=0.2]{<->}(5,0)(0,0)(0,5)
-\psline[arrowinset=0.05,arrowsize=0.1]{<->}(5,0)(0,0)(0,5)
-\uput[u](0,5){$y$}\uput[r](5,0){$x$}
-\end{pspicture}
- \end{center}
-Il s'obtient très facilement à partir de l'étude précédente si on sait que dans ce cas :
-\[
-v_0=\sqrt{\frac{\mathcal{G}M}{r_0}}
-\]
-\begin{verbatim}
-\pstVerb{
-    /GM 1 def
-    /theta0 60 def
-    /r0 5 def
-    /x0 r0 theta0 cos mul def
-    /y0 r0 theta0 sin mul def
-    /v0 GM r0 div sqrt def
-    /v0x v0 theta0 sin mul neg def
-    /v0y v0 theta0 cos mul def }%
-\end{verbatim}
-
-\end{document}
+\documentclass{article}\r
+\usepackage[a4paper,margin=2cm]{geometry}\r
+\usepackage[latin1]{inputenc}\r
+\usepackage{pstricks-add,pst-eqdf,pst-func}\r
+\usepackage{array,amsmath}\r
+\usepackage{listings}\r
+\newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.125,linecolor={[rgb]{1 0.5 0}}}\r
+\newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}\r
+\newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0.5 0}}}\r
+%%%%%%%%%%%%%%%%%%\r
+\title{Interaction gravitationnelle avec PSTricks}\r
+\date{19 juin 2\,012}\r
+\begin{document}\r
+\maketitle\r
+\def\eqsatellite{%\r
+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)}%\r
+\section{Mise en orbite d'un satellite}\r
+\begin{figure}[htbp]\r
+  \begin{center}\r
+\psset{unit=2}\r
+\begin{pspicture}(-6,-2)(2,5)\r
+\pstVerb{\r
+    /GM 1 def\r
+    /theta0 -35 def\r
+    /r0 1 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 1.3 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /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\r
+%%%%%%%%%%%%%%\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+    /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def\r
+% vitesse à l'apogée\r
+    /vA GM par div sqrt 1 exc sub mul def\r
+% vitesse au périgée\r
+    /vP GM par div sqrt 1 exc add mul def\r
+% coordonnées de vA\r
+    /vAx vA theta0 90 add cos mul neg def\r
+    /vAy vA theta0 90 add sin mul neg def\r
+% coordonnées de vP\r
+    /vPx vP theta0 90 add cos mul def\r
+    /vPy vP theta0 90 add sin mul def\r
+% Apogée\r
+    /xA par 1 exc sub div theta0 cos mul neg def\r
+    /yA par 1 exc sub div theta0 sin mul neg def\r
+% Périgée\r
+   /xP par 1 exc add div theta0 cos mul def\r
+   /yP par 1 exc add div theta0 sin mul def\r
+% Centre\r
+   /xO xP xA add 2 div def\r
+   /yO yP yA add 2 div def\r
+}%\r
+\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}\r
+\psdot[dotstyle=+](0,0)\r
+\uput[d](0,0){O}\r
+\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-2)(2,5)\r
+\rput(1,4){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{0000}s}}\r
+\parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def\r
+ radius t cos mul\r
+ radius t sin mul}\r
+\pscircle*(!x0 y0){0.05}\r
+\pnode(!xP yP){P} % périgée\r
+\pnode(!xA yA){A} % Apogée\r
+\pnode(!xO yO){O} % centre\r
+\rput(!x0 y0){\psline[style=vecteurB]{->}(!v0x v0y)\uput[ur](!v0x v0y){$\overrightarrow{v_0}$}}\r
+\psline[linestyle=dashed](A)(P)\r
+\pscircle*(A){0.05}\r
+\uput[ul](A){$A$}\r
+\uput[dr](P){$P$}\r
+% position du satellite à un instant quelconque\r
+\pstVerb{/theta_i 170 def\r
+ /radius par 1 exc theta_i theta0 sub cos mul add div def\r
+ /xS radius theta_i cos mul def\r
+ /yS radius theta_i sin mul def\r
+ /ux theta_i cos 1 mul def\r
+ /uy theta_i sin 1 mul def\r
+ /xi xS ux sub def\r
+ /yi yS uy sub def\r
+ /xi2 xS ux 2 div sub def\r
+ /yi2 yS uy 2 div sub def}%\r
+\pnode(!xi2 yi2){Mi2}\r
+\pnode(!xi yi){Mi}\r
+\pnode(!xS yS){S}\r
+\pscircle*(S){0.05}\r
+\psline[linestyle=dotted](S)(0,0)\r
+\psline[style=vecteurA]{->}(S)(Mi)\r
+\uput[l](S){S}\r
+\uput[u](Mi2){$\overrightarrow{F}$}\r
+\psarcn{->}(0,0){0.4}{0}{!theta0}\r
+\uput{0.5}[!theta0 2 div](0,0){$\theta_0$}\r
+\rput(A){\psline[style=vecteurB]{->}(!vAx vAy)}\r
+%\rput(P){\psline[style=vecteurB]{->}(!vPx vPy)}\r
+\psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,5)\r
+\uput[r](0,4.9){$y$}\uput[u](1.9,0){$x$}\r
+\psdot(O)\uput[r](O){$\Omega$}\r
+\rput{!90 theta0 add}(O){\psline[linestyle=dashed](!b_2 neg 0)(!b_2 0)}\r
+\end{pspicture}\r
+\caption{Mouvement d'un satellite}\r
+  \end{center}\r
+\end{figure}\r
+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 :\r
+\[\r
+\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}\r
+\]\r
+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 :\r
+\[\r
+r=\frac{p}{1+\mathrm{e}\cos(\theta-\theta_0)}\r
+\]\r
+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.\r
+\[\r
+\mathrm{e}=\sqrt{1+\frac{2\mathcal{E}L^2}{mK^2}}\qquad p=\frac{L^2}{mK}\r
+\]\r
+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 :\r
+\[\r
+L=mr_0v_0\qquad \mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2\r
+\]\r
+En remplaçant $L$ et $\mathcal{E}$, on obtient pour l'excentricité et le paramètre les expressions suivantes :\r
+\[\r
+\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)}\r
+\]\r
+\[\r
+p=\frac{v_0^2r_0^2}{\mathcal{G}M}\r
+\]\r
+On se limite au cas du mouvement elliptique, avec, en conséquence, la condition :\r
+\[\r
+\mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2 < 0\r
+\]\r
+Le demi-grand axe $a$, le demi-petit axe $b$ sont :\r
+\[\r
+a=\frac{p}{1-\mathrm{e}^2}\qquad b=\frac{p}{\sqrt{1-\mathrm{e}^2}}\r
+\]\r
+La période $T$ qui obéit à la troisième loi de Képler :\r
+\[\r
+T^2=\frac{4\pi^2a^3}{\mathcal{G}M}\r
+\]\r
+La vitesse, en un point de l'ellipse, se calcule par :\r
+\[\r
+v^2=\mathcal{G}M\Big( \frac{2}{r}-\frac{1}{a}\Big)\r
+\]\r
+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 :\r
+\[\r
+v_P=\sqrt{\frac{\mathcal{G}M}{p}}(1+\mathrm{e})\qquad v_A=\sqrt{\frac{\mathcal{G}M}{p}}(1-\mathrm{e})\r
+\]\r
+\r
+\section{L'étude avec PSTricks}\r
+\subsection{La trajectoire}\r
+\def\parametres{\r
+    /GM 1 def % 4e14 def % GxM\r
+    /x0 6.5e6 def % position initiale\r
+    /y0 0 def\r
+    /vx0 0 def    % vitesse initiale\r
+    /vy0 1e4 def\r
+}\r
+%  x0   y0   x'0   y'0\r
+% y[0] y[1] y[2]  y[3]\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+\ddot{x}&=&-\dfrac{GM}{r^3}x\\[1em]\r
+\ddot{y}&=&-\dfrac{GM}{r^3}y\\\r
+\end{array}\r
+\right.\r
+\label{eq1}\r
+\qquad r=\sqrt{x^2+y^2}\r
+\]\r
+\r
+On peut dessiner la trajectoire du satellite et de ses caractéristiques de deux façons :\r
+\begin{itemize}\r
+  \item par l'utilisation de \verb+\parametricplot+ ;\r
+  \item ou celle de \verb+\psplotDiffEqn+.\r
+\end{itemize}\r
+\verb+\parametricplot+ utilise l'expression exacte de l'équation de la trajectoire en coordonnées polaires :\r
+\begin{verbatim}\r
+   \parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{%\r
+     /radius par 1 exc t theta0 sub cos mul add div def\r
+      radius t cos mul\r
+      radius t sin mul}\r
+\end{verbatim}\r
+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}.\r
+\begin{verbatim}\r
+\pstVerb{\r
+    /GM 1 def\r
+    /theta0 -45 def\r
+    /r0 0.5 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 1.92 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul\r
+         v0 dup mul mul sub GM dup mul div 2 mul add sqrt def\r
+%demi-grand axe\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+%demi-petit axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+% période\r
+    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def\r
+}%\r
+\end{verbatim}\r
+\verb+\psplotDiffEqn+ utilise les équations différentielles du mouvement, en notation algébrique :\r
+\begin{verbatim}\r
+%  x0   y0   x'0   y'0\r
+% y[0] y[1] y[2]  y[3]\r
+\def\eqsatellite{%\r
+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)}\r
+ \psplotDiffEqn[unit=2,whichabs=0,whichord=1,%\r
+                linecolor=blue,linewidth=0.1,%\r
+                method=rk4,plotpoints=1000,%\r
+                algebraic]{0}{37.8}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\end{verbatim}\r
+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.\r
+%\r
+\begin{center}\r
+\begin{pspicture}(-12,-2)(4,10)\r
+\pstVerb{\r
+    /GM 1 def\r
+    /theta0 -45 def\r
+    /r0 0.5 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 1.92 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /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\r
+%%%%%%%%%%%%%%\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def\r
+}%\r
+\psframe*[linecolor=white](-3,-0.2)(0,0.2)\r
+\pscircle[fillcolor=blue!50,fillstyle=solid](0,0){0.5}\r
+\psgrid[unit=2,subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-1)(2,5)\r
+\rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}\r
+  \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}%\r
+\parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def\r
+ radius t cos mul\r
+ radius t sin mul}\r
+\psdot[unit=2,dotsize=0.12](!x0 y0)\r
+\rput(!x0 2 mul y0 2 mul){\psline[style=vecteurC]{->}(!v0x v0y)}\r
+\end{pspicture}\r
+\end{center}\r
+\subsection{La vitesse}\r
+\verb+\psplotDiffEqn+ permet de voir comment varie la vitesse sur l'ellipse :\r
+\begin{verbatim}\r
+%  x0   y0   x'0   y'0\r
+% y[0] y[1] y[2]  y[3]\r
+ \psplotDiffEqn[xunit=0.2,yunit=5,%\r
+                plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,\r
+                linecolor=red,method=rk4,plotpoints=1000,\r
+                algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\end{verbatim}\r
+\begin{center}\r
+\begin{pspicture}(0,-1)(10,7)\r
+\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]\r
+\pstVerb{/GM 1 def\r
+    /theta0 -35 def\r
+    /r0 1 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 1.3 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def}%\r
+\psplotDiffEqn[xunit=0.2,yunit=5,\r
+  plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,\r
+  linecolor=red,method=rk4,plotpoints=1000,algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%\r
+ \multido{\i=1+1,\I=5+5}{9}{\uput[u](\i,0){\I}}\r
+\pnode(! 36.4 5 div 0){P}\r
+\psdot(P)\uput[d](P){Périgée}\r
+\psline[linestyle=dashed](P)(! 36.4 5 div 7)\r
+\pnode(! 36.4 10 div 0){A}\r
+\psdot(A)\uput[d](A){Apogée}\r
+\psline[linestyle=dashed](A)(! 36.4 10 div 7)\r
+\uput[l](0,6.5){$v$}\r
+\psline[arrowinset=0.1,arrowsize=0.2]{<->}(10,0)(0,0)(0,7)\r
+\uput[u](10,0){$t$(s)}\r
+\uput[ul](0,0){0}\r
+\end{pspicture}\r
+\end{center}\r
+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 :\r
+\[\r
+\overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}\r
+\]\r
+$\dot{\theta}$ et $\dot{r}$ s'obtiennent par les relations suivantes :\r
+\[\r
+\dot{\theta}=\frac{r_0v_0}{r^2}\r
+\]\r
+\[\r
+\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)\r
+\]\r
+La chaîne de calculs est la suivante : $\theta\Longrightarrow r\Longrightarrow \dot{\theta}\Longrightarrow \dot{r}\Longrightarrow \overrightarrow{v}$.\r
+\r
+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.\r
+\r
+On sauve successivement le tableau des positions et celui des vitesses.\r
+\begin{verbatim}\r
+\psequadiff[method=rk4,plotpoints=1000,\r
+            algebraic,\r
+            whichabs=0,whichord=1,\r
+            tabname=XiYi]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\end{verbatim}\r
+\begin{verbatim}\r
+\psequadiff[method=rk4,plotpoints=1000,\r
+            algebraic,\r
+            whichabs=2,whichord=3,\r
+            tabname=vxvy]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\end{verbatim}\r
+\r
+Pour ensuite dessiner la trajectoire et les vecteurs-vitesse.\r
+\begin{verbatim}\r
+%\listplot[unit=1]{vxvy aload pop}\r
+% on dessine la vitesse un point sur 100\r
+\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}\r
+\multido{\i=0+100}{20}{%\r
+\pstVerb{/vX vxvy \i\space get def\r
+         /vY vxvy \i\space 1 add get def\r
+         /xi XiYi \i\space get def\r
+         /yi XiYi \i\space 1 add get def}%\r
+\rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}\r
+\end{verbatim}\r
+\begin{center}\r
+\begin{pspicture}(-10,-10)(6,7)\r
+\psset{unit=2}%\r
+\pstVerb{/GM 1 def\r
+    /theta0 30 def\r
+    /r0 2 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 0.85 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /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\r
+%%%%%%%%%%%%%%\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def}%\r
+\rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}\r
+\psequadiff[method=rk4,\r
+            plotpoints=1000,\r
+            algebraic,\r
+            whichabs=0,\r
+            whichord=1,\r
+            tabname=XiYi\r
+%           ,saveData,filename=XiYi.dat\r
+]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\listplot{XiYi aload pop}\r
+\psequadiff[method=rk4,\r
+            plotpoints=1000,\r
+            algebraic,\r
+            whichabs=2,\r
+            whichord=3,\r
+            tabname=vxvy\r
+%           ,saveData,filename=vxvy.dat\r
+]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%\r
+%\listplot[unit=1]{vxvy aload pop}\r
+% on dessine la vitesse un point sur 100\r
+\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}\r
+\multido{\i=0+100}{20}{%\r
+\pstVerb{/vX vxvy \i\space get def\r
+         /vY vxvy \i\space 1 add get def\r
+         /xi XiYi \i\space get def\r
+         /yi XiYi \i\space 1 add get def}%\r
+\rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}\r
+\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(3,3)\r
+\end{pspicture}\r
+\end{center}\r
+%\r
+\section{Mouvement circulaire}\r
+ \begin{center}\r
+\psset{unit=1}\r
+\begin{pspicture}(-5,-5.5)(5,5.5)\r
+\pstVerb{\r
+    /GM 1 def\r
+    /theta0 60 def\r
+    /r0 5 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 GM r0 div sqrt def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /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\r
+%%%%%%%%%%%%%%\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+    /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def\r
+% vitesse à l'apogée\r
+    /vA GM par div sqrt 1 exc sub mul def\r
+% vitesse au périgée\r
+    /vP GM par div sqrt 1 exc add mul def\r
+% coordonnées de vA\r
+    /vAx vA theta0 90 add cos mul neg def\r
+    /vAy vA theta0 90 add sin mul neg def\r
+% coordonnées de vP\r
+    /vPx vP theta0 90 add cos mul def\r
+    /vPy vP theta0 90 add sin mul def\r
+}%\r
+\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.75}\r
+\psdot[dotstyle=+](0,0)\r
+\uput[d](0,0){O}\r
+\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(5,5)\r
+\rput(0,-2){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{000000}s}}\r
+\parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def\r
+ radius t cos mul\r
+ radius t sin mul}\r
+\pscircle*(!x0 y0){0.1}\r
+%\pnode(!par 1 exc add div theta0 cos mul par 1 exc add div theta0 sin mul){P} % périgée\r
+%\pnode(!par 1 exc sub div theta0 cos mul neg par 1 exc sub div theta0 sin mul neg){A} % Apogée\r
+\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}$}}\r
+\uput[ur](!x0 y0){$M_0$}\r
+\psline[linestyle=dashed](0,0)(!x0 y0)\r
+% position du satellite à un instant quelconque\r
+\pstVerb{/theta_i 170 def\r
+ /radius par 1 exc theta_i theta0 sub cos mul add div def\r
+ /xS radius theta_i cos mul def\r
+ /yS radius theta_i sin mul def\r
+ /ux theta_i cos 1 mul def\r
+ /uy theta_i sin 1 mul def\r
+ /xi xS ux sub def\r
+ /yi yS uy sub def\r
+ /xi2 xS ux 2 div sub def\r
+ /yi2 yS uy 2 div sub def}%\r
+\pnode(!xi2 yi2){Mi2}\r
+\pnode(!xi yi){Mi}\r
+\pnode(!xS yS){S}\r
+\pscircle*(S){0.1}\r
+\psline[linestyle=dotted](S)(0,0)\r
+\psline[style=vecteurC]{->}(S)(Mi)\r
+\uput[l](S){S}\r
+\uput[u](Mi2){$\overrightarrow{F}$}\r
+\psarc{->}(0,0){1}{0}{!theta0}\r
+\uput{1.1}[!theta0 2 div](0,0){$\theta_0$}\r
+\psline[arrowinset=0.1,arrowsize=0.2]{<->}(5,0)(0,0)(0,5)\r
+\psline[arrowinset=0.05,arrowsize=0.1]{<->}(5,0)(0,0)(0,5)\r
+\uput[u](0,5){$y$}\uput[r](5,0){$x$}\r
+\end{pspicture}\r
+ \end{center}\r
+Il s'obtient très facilement à partir de l'étude précédente si on sait que dans ce cas :\r
+\[\r
+v_0=\sqrt{\frac{\mathcal{G}M}{r_0}}\r
+\]\r
+\begin{verbatim}\r
+\pstVerb{\r
+    /GM 1 def\r
+    /theta0 60 def\r
+    /r0 5 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 GM r0 div sqrt def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def }%\r
+\end{verbatim}\r
+\section{L'hodographe du mouvement du satellite}\r
+On rappelle qu'en coordonnées polaires le vecteur-vitesse s'écrit :\r
+\[\r
+\overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}\r
+\]\r
+\[\r
+\dot{\theta}=\frac{r_0v_0}{r^2}\r
+\]\r
+\[\r
+\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)\r
+\]\r
+On peut exprimer $v$ uniquement en fonction de $\theta$ :\r
+\[\r
+\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}\r
+\]\r
+Ses composantes dans la base $(\overrightarrow{u_r},u_\theta)$ sont :\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+v_r&=&\dfrac{r_0v_0}{p}\mathrm{e}\sin(\theta-\theta_0)\\[1em]\r
+v_\theta&=&\dfrac{r_0v_0}{p}\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)\r
+\end{array}\r
+\right.\r
+\label{eqv1}\r
+\]\r
+On repasse aux coordonnées cartésiennes par une rotation d'angle $(-\theta)$.\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+\dot{x}&=&v_r\cos\theta-v_\theta\sin\theta\\[1em]\r
+\dot{y}&=&v_r\sin\theta+v_\theta\cos\theta\r
+\end{array}\r
+\right.\r
+\label{eqv2}\r
+\]\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+\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]\r
+\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\r
+\end{array}\r
+\right.\r
+\label{eqv3}\r
+\]\r
+En utilisant les relations trigonométriques de soustraction :\r
+\[\r
+\sin\alpha\cos\beta-\cos\alpha\sin\beta=\sin(\alpha-\beta)\qquad \cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta)\r
+\]\r
+On obtient :\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+\dfrac{p}{r_0v_0}\dot{x}&=&-\mathrm{e}\sin\theta_0-\sin\theta\\[1em]\r
+\dfrac{p}{r_0v_0}\dot{y}&=&\hphantom{-}\mathrm{e}\cos\theta_0+\cos\theta\r
+\end{array}\r
+\right.\r
+\label{eqv4}\r
+\]\r
+\[\r
+\left\{\r
+\begin{array}{rcl}\r
+\dot{x}+\dfrac{r_0v_0}{p}\mathrm{e}\sin\theta_0&=&-\dfrac{r_0v_0}{p}\sin\theta\\[1em]\r
+\dot{y}-\dfrac{r_0v_0}{p}\mathrm{e}\cos\theta_0&=&\dfrac{r_0v_0}{p}\cos\theta\r
+\end{array}\r
+\right.\r
+\label{eqv5}\r
+\]\r
+L'équation de l'hodographe :\r
+\[\r
+\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\r
+\]\r
+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.\r
+\begin{center}\r
+\begin{pspicture}(-12,-2)(4,8)\r
+\psset{unit=2}%\r
+\uput[r](0,4){$y$}\uput[u](2,0){$x$}\r
+\pstVerb{/GM 1 def\r
+    /theta0 -35 def\r
+    /r0 1 def\r
+    /x0 r0 theta0 cos mul def\r
+    /y0 r0 theta0 sin mul def\r
+    /v0 1.3 def\r
+    /v0x v0 theta0 sin mul neg def\r
+    /v0y v0 theta0 cos mul def\r
+    /Lc r0 v0 mul def % moment cinetique\r
+    /par Lc dup mul GM div def % paramètre de l'ellipse\r
+% excentricité\r
+    /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\r
+%%%%%%%%%%%%%%\r
+    /a_2 par 1 exc dup mul sub div def % demi-grand axe\r
+    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe\r
+    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def\r
+% rayon de l'hodographe\r
+   /rH r0 v0 mul par div def\r
+% les coordonnées du centre\r
+   /xCH rH exc mul theta0 sin mul neg def\r
+   /yCH rH exc mul theta0 cos mul def}%\r
+\psequadiff[method=rk4,\r
+            plotpoints=1000,\r
+            algebraic,\r
+            whichabs=0,\r
+            whichord=1,\r
+            tabname=XiYi,\r
+%            saveData,filename=XiYi.dat\r
+]{0}{36.5}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\listplot{XiYi aload pop}\r
+\psequadiff[method=rk4,\r
+            plotpoints=1000,\r
+            algebraic,\r
+            whichabs=2,\r
+            whichord=3,\r
+            tabname=vxvy,\r
+%            saveData,filename=vxvy.dat\r
+]{0}{36.5}{x0 y0 v0x v0y}{\eqsatellite}%\r
+\listplot[unit=1,linecolor=red,linewidth=0.075]{vxvy aload pop}\r
+% on dessine la vitesse un point sur 50\r
+\multido{\i=0+50}{40}{%\r
+\pstVerb{ /vX vxvy \i\space get def\r
+          /vY vxvy  \i\space 1 add get def\r
+          /xi XiYi \i\space get def\r
+          /yi XiYi \i\space 1 add get def}%\r
+\rput(!xi yi){\psline[unit=1,linecolor=blue]{->}(!vX vY)}}\r
+\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.2}\r
+\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-2)(2,4)\r
+\psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,4)\r
+\pscircle(!xCH yCH){!rH}\r
+\psdot[dotstyle=+](!xCH yCH)\r
+\end{pspicture}\r
+\end{center}\r
+\end{document}\r

Licence Creative Commons Les fichiers de Syracuse sont mis à disposition (sauf mention contraire) selon les termes de la
Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Partage dans les Mêmes Conditions 4.0 International.