\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage[latin1]{inputenc}
\usepackage{pst-eqdf,pst-func}
\usepackage{array,amsmath}
\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{L'hodographe du mouvement d'un satellite}
\date{23 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)}%
On rappelle qu'en coordonnées polaires le vecteur-vitesse s'écrit :
\[
\overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}
\]
\[
\dot{\theta}=\frac{r_0v_0}{r^2}
\]
\[
\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)
\]
On peut exprimer $v$ uniquement en fonction de $\theta$ :
\[
\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}
\]
Ses composantes dans la base $(\overrightarrow{u_r},u_\theta)$ sont :
\[
\left\{
\begin{array}{rcl}
v_r&=&\dfrac{r_0v_0}{p}\mathrm{e}\sin(\theta-\theta_0)\\[1em]
v_\theta&=&\dfrac{r_0v_0}{p}\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)
\end{array}
\right.
\label{eqv1}
\]
On repasse aux coordonnées cartésiennes par une rotation d'angle $(-\theta)$.
\[
\left\{
\begin{array}{rcl}
\dot{x}&=&v_r\cos\theta-v_\theta\sin\theta\\[1em]
\dot{y}&=&v_r\sin\theta+v_\theta\cos\theta
\end{array}
\right.
\label{eqv2}
\]
\[
\left\{
\begin{array}{rcl}
\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]
\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
\end{array}
\right.
\label{eqv3}
\]
En utilisant les relations trigonométriques de soustraction :
\[
\sin\alpha\cos\beta-\cos\alpha\sin\beta=\sin(\alpha-\beta)\qquad \cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta)
\]
On obtient :
\[
\left\{
\begin{array}{rcl}
\dfrac{p}{r_0v_0}\dot{x}&=&-\mathrm{e}\sin\theta_0-\sin\theta\\[1em]
\dfrac{p}{r_0v_0}\dot{y}&=&\hphantom{-}\mathrm{e}\cos\theta_0+\cos\theta
\end{array}
\right.
\label{eqv4}
\]
\[
\left\{
\begin{array}{rcl}
\dot{x}+\dfrac{r_0v_0}{p}\mathrm{e}\sin\theta_0&=&-\dfrac{r_0v_0}{p}\sin\theta\\[1em]
\dot{y}-\dfrac{r_0v_0}{p}\mathrm{e}\cos\theta_0&=&\dfrac{r_0v_0}{p}\cos\theta
\end{array}
\right.
\label{eqv5}
\]
L'équation de l'hodographe :
\[
\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
\]
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.
\begin{center}
\begin{pspicture}(-10,-2)(4,8)
\psset{unit=2}%
\uput[r](0,4){$y$}\uput[u](2,0){$x$}
\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 dup mul a_2 3 exp mul GM div sqrt mul def
% rayon de l'hodographe
   /rH r0 v0 mul par div def
% les coordonnées du centre
   /xCH rH exc mul theta0 sin mul neg def
   /yCH rH exc mul theta0 cos mul def}%
\psequadiff[method=rk4,
            plotpoints=1000,
            algebraic,
            whichabs=0,
            whichord=1,
            tabname=XiYi,
            saveData,filename=XiYi.dat]{0}{36.5}{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}{36.5}{x0 y0 v0x v0y}{\eqsatellite}%
\listplot[unit=1,linecolor=red,linewidth=0.075]{vxvy aload pop}
% on dessine la vitesse un point sur 50
\multido{\i=0+50}{40}{%
\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[unit=1,linecolor=blue]{->}(!vX vY)}}
\pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.2}
\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-2)(2,4)
\psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,4)
\pscircle(!xCH yCH){!rH}
\psdot[dotstyle=+](!xCH yCH)
\end{pspicture}
\end{center}
\newpage
\begin{center}
\begin{pspicture}(-10,-10)(6,6)
\psset{unit=2}%
\uput[r](0,3){$y$}\uput[u](3,0){$x$}
\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=2,linecolor=red]{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)
\psline[arrowinset=0.05,arrowsize=0.1]{<->}(3,0)(0,0)(0,3)
\end{pspicture}
\end{center}
\end{document}
