\documentclass[french, a4paper, 12pt, reqno]{amsart}% article ams
\usepackage[french]{babel}% prise en compte les particularites de la typographie francaise. Remplace francais et frenchb
\usepackage[T1]{fontenc}% codage de caractères T1, nouvelle norme LaTeX. Evite l'erreur [OT1] mais passe en Bitmap
\usepackage{lmodern} % fonte latin modern. Re-passe de Bitmap à vectoriel
\usepackage{amssymb, latexsym}% symboles mathematiques
\usepackage{pst-all}% dessiner*
\usepackage{xspace}% gestion des espaces et demi-espaces
\usepackage[left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}% mise en page
\usepackage{float}% begin figure[H] (ne marche pas avec subfigure, conserver [h!])
\usepackage{empheq}% mise en évidence des formules importantes
\definecolor{mygrey}{rgb}{.9, .9, .9}
\newcommand*\greybox[1]{\colorbox{mygrey}{\hspace{1em}#1\hspace{1em}}}
\newtheorem*{remarque}{Remarque}
\newenvironment{rmq}{\begin{remarque}\small}{\end{remarque}}
\usepackage{fancyhdr}% numérotation des pages
\fancyhf{}
\pagestyle{fancy}
\fancyhead[LE,RO]{\bfseries\thepage}% position des numeros de pages
\fancyfoot{}% rien en pied de page
\fancypagestyle{plain}{\renewcommand{\headrulewidth}{0pt}\fancyhf{}}% supprime le numéro de page des nouveaux chapitres
\renewcommand{\headrulewidth}{0pt}% pas de filet en haut des pages
\setlength{\headheight}{16pt}
\usepackage{hyperref}
\hypersetup{
    pdftoolbar=true,        % affiche la barre d'outils d'Acrobat
    pdfmenubar=true,        % affiche le menu d'Acrobat
    pdftitle={Le problème de Kepler},    % titre
    pdfauthor={Olivier Castéra},     % auteur
    pdfstartview= FitH,		%La page prend toute la largeur
    colorlinks=true,       % false: liens encadres ; true: liens colores
    linkcolor=blue,        % couleur des liens internes
    }

\title{Le problème de Kepler}
\author{Olivier Castéra}
\email{o.castera@free.fr}
\urladdr{http://o.castera.free.fr/}
\date{\today}
\usepackage{titlesec}% formatage des titres
\titleformat{\section}[block]{\vspace{\baselineskip}\large\scshape}{\thesection}{.7cm}{}[\hrule\vspace{.75\baselineskip}]
\titleformat{\subsection}[block]{\vspace{.75\baselineskip}\bfseries}{\thesubsection}{.5cm}{}[]
\titleformat{\subsubsection}[block]{\vspace{.75\baselineskip}\itshape}{\thesubsubsection}{.2cm}{}[]
\renewcommand{\numberline}[1]{#1~~}% Problème de numérotation de la ToC

\newtheorem{theoreme}{Théorème}[section]
\newtheorem{definition}{Définition}[section]

\def\grad{\mathop{\boldsymbol{\rm grad}}\nolimits}

\begin{document}
\begin {abstract}
Il s'agit de déterminer l'équation de l'orbite d'une pla\-nè\-te à partir de la mécanique de Newton,
en ne considérant que la seule force de gravitation, en négligeant l'interaction gravitationnelle des planètes entre elles,
et en supposant les astres ponctuels.
Ce problème est aussi connu sous le nom de \og problème des deux corps\fg.
\end {abstract}
\maketitle% affiche le titre
\thispagestyle{empty}% supprime l'affichage du numéro de la première page
\tableofcontents% affiche la table des matières (ToC)
\section{Coordonnées polaires}
Etant donné que l'on traite ici de la révolution d'une planète autour du Soleil,
il est avantageux de passer en coordonnées cylindriques $(\rho,\theta,z)$.
Les axes des $z$ en coordonnées rectangulaires $(x,y,z)$ et en cylindriques sont supposés confondus.
Nous nous intéressons ici aux coordonnées polaires $(\rho,\theta)$ prises dans le plan $(x,y)$.
\subsection{Expression des vecteurs de base de la base polaire}

\begin{figure}[H]
\begin{pspicture}(-4,-1)(6,4)%aire reservee a l'image
%\psgrid(-4,-1)(6,4)
\psline{->}(0,0)(5,0)
\psline{->}(0,0)(0,3)
\rput(5.3,0){$x$}
\rput(0,3.3){$y$}
\rput(-.3,-.3){$O$}
\rput(.8,-.3){$\boldsymbol i$}
\rput(-.3,.8){$\boldsymbol j$}
\rput(3,.7){$M$}
\rput(1.8,1){$\rho$}
\SpecialCoor
\psdot(3;20)
\psline{-}(0,0)(3;20)
\psline{->}(0,0)(1,0)
\psline{->}(0,0)(0,1)
\rput(3;20){
\psline[linecolor=red]{->}(0;0)(1;20)
\psline[linecolor=red]{->}(0;0)(1;110)
\rput(1.3;110){\red $\boldsymbol e_{\theta}$}
\rput(1.3;20){\red $\boldsymbol e_{\rho}$}
}
\psarc{<->}(0;0){1.5}{0}{20}
\rput(1.8,.3){$\theta$}
\psarc{->}(0;0){4.8}{5}{15}
\rput(5.1;10){+}
\end{pspicture}
\caption{Vecteurs de la base polaire}
\label{bp}
\end{figure}

\subsubsection*{Première méthode}
En se servant de la figure~\ref{bp}, nous pouvons exprimer les vec\-teurs de la base polaire en fonction de ceux de la base cartésienne :
\begin{equation*}
\left\{
\begin{aligned}
&\boldsymbol e_\rho=\cos\theta \boldsymbol i+\sin\theta \boldsymbol j\\
&\boldsymbol e_\theta=-\sin\theta \boldsymbol i+\cos\theta \boldsymbol j
\end{aligned}
\right.
\end{equation*}
\subsubsection*{Deuxième méthode}
Soit $M$ un point de coordonnées polaires $(\rho,\theta)$ et de coordon\-nées cartésiennes $(x,y)$.
Le changement de coordonnées de polaires vers cartésiennes s'écrit~:
\begin{equation*}
\left\{
\begin{aligned}
&x=\rho\cos\theta\\
&y=\rho\sin\theta
\end{aligned}
\right.
\end{equation*}
Déterminons les vecteurs de base $\boldsymbol e_\rho$ et $\boldsymbol e_\theta $
en différentiant le rayon vecteur $\boldsymbol r=\boldsymbol{OM}$
exprimé dans la base cartésienne $(\boldsymbol i,\boldsymbol j)$ en fonction des coordonnées polaires $(\rho,\theta)$~:
\begin{align*}
\boldsymbol r(x,y)&=x\boldsymbol i+y\boldsymbol j\\
\boldsymbol r(\rho,\theta)&=\rho\cos\theta\,\boldsymbol i+\rho\sin\theta\,\boldsymbol j\\
d\boldsymbol r(\rho,\theta)&=\left(\frac{\partial\boldsymbol r}{\partial\rho}\right)_\theta d\rho
+\left(\frac{\partial\boldsymbol r}{\partial\theta}\right)_\rho d\theta\\
&=(\cos\theta\,\boldsymbol i+\sin\theta\,\boldsymbol j)d\rho+\rho(-\sin\theta\,\boldsymbol i+\cos\theta\,\boldsymbol j)d\theta
\end{align*}
Les vecteurs unitaires de la base polaire ont alors pour expression,
\begin{equation*}
\left\{
\begin{aligned}
&\boldsymbol e_\rho=\cos\theta \boldsymbol i+\sin\theta \boldsymbol j\\
&\boldsymbol e_\theta=-\sin\theta \boldsymbol i+\cos\theta \boldsymbol j
\end{aligned}
\right.
\end{equation*}
et l'on a~:
\begin{align*}
d\boldsymbol r(\rho,\theta)=\boldsymbol e_\rho\,d\rho+\rho\boldsymbol e_\theta \,d\theta
\end{align*}
\subsection{Expression du vecteur position}
\begin{align*}
\boldsymbol r(\rho,\theta)&=\rho\cos\theta\,\boldsymbol i+\rho\sin\theta\,\boldsymbol j\\
&=\rho(\cos\theta\,\boldsymbol i+\sin\theta\,\boldsymbol j)\\
&=\rho \boldsymbol e_\rho
\end{align*}
Dans la base polaire $(\boldsymbol e_\rho,\boldsymbol e_\theta )$ et en coordonnées polaires $(\rho,\theta)$, le vecteur position s'écrit~:
\begin{empheq}[box=\greybox]{align*}
\boldsymbol r=\rho \boldsymbol e_\rho
\end{empheq}
\subsection{Dérivée des vecteurs de base}
Nous avons besoin de la dérivée des vecteurs de base pour exprimer la vitesse et l'accélération~:
\begin{equation*}
\left\{
\begin{aligned}
&\frac{d\boldsymbol e_\rho}{dt}=-\sin\theta\,\frac{d\theta}{dt}\,\boldsymbol i+\cos\theta\,\frac{d\theta}{dt}\,\boldsymbol j\\
&\frac{d\boldsymbol e_\theta}{dt}=-\cos\theta\,\frac{d\theta}{dt}\,\boldsymbol i-\sin\theta\,\frac{d\theta}{dt}\,\boldsymbol j
\end{aligned}
\right.
\qquad\Rightarrow\qquad
\left\{
\begin{aligned}
&\dot{\boldsymbol e}_\rho=\dot\theta\boldsymbol e_\theta\\
&\dot{\boldsymbol e}_\theta=-\dot\theta\boldsymbol e_\rho
\end{aligned}
\right.
\end{equation*}
\subsection{Expression du vecteur vitesse}
Le vecteur vitesse est la dérivée première du vecteur position par rapport au temps~:
\begin{align*}
\boldsymbol v&=\frac{d\boldsymbol r}{dt}\\
&=\frac{d}{dt}\left(\rho\boldsymbol e_\rho \right)\\
&=\dot\rho\boldsymbol e_\rho+\rho\dot{\boldsymbol e}_\rho
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{vit_pol}
\boldsymbol v&=\dot\rho\boldsymbol e_\rho+\rho\dot\theta\boldsymbol e_\theta 
\end{empheq}
\subsection{Expression du vecteur accélération}
Le vecteur accélération est la dérivée première du vecteur vitesse par rapport au temps~:
\begin{align*}
\boldsymbol{a}&=\frac{d\boldsymbol v}{dt}\\
&=\frac{d}{dt}\left(\dot\rho\boldsymbol e_\rho+\rho\dot\theta\boldsymbol e_\theta\right)\\
&=\ddot{\rho}\boldsymbol e_\rho+\dot\rho\dot{\boldsymbol e}_\rho+\dot\rho\dot\theta\boldsymbol e_\theta
+\rho\ddot\theta\boldsymbol e_\theta +\rho\dot\theta\dot{\boldsymbol e}_\theta\\
&=\ddot{\rho}\boldsymbol e_\rho+\dot\rho\dot\theta\boldsymbol e_\theta+\dot\rho\dot\theta\boldsymbol e_\theta
+\rho\ddot\theta\boldsymbol e_\theta-\rho\dot\theta^2\boldsymbol e_\rho
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{acc_pol}
\boldsymbol{a}
&=\left(\ddot{\rho}-\rho\dot\theta^2\right)\boldsymbol e_\rho+\left(2\dot\rho\dot\theta+\rho\ddot\theta\right)\boldsymbol e_\theta 
\end{empheq}
\section{Relation fondamentale de la dynamique}
\subsection{Centre d'inertie}
Soient $m_1$ et $m_2$ les masses des astres.
Soit $O$ un point quelconque fixe dans un référentiel galiléen $\mathcal R$,
notons $\boldsymbol{OA_1}$ et $\boldsymbol{OA_2}$ les vecteurs position des astres dans $\mathcal R$.
\begin{definition}
Le point $G$ fixe dans $\mathcal R$ tel que
\begin{align*}
&m_1\,\boldsymbol{GA_1}+m_2\,\boldsymbol{GA_2}=\boldsymbol 0
\end{align*}
est appelé \textit{centre de masse} ou \textit{centre d'inertie} du système.
\end{definition}
\begin{definition}
Le référentiel de centre $G$ est appelé référentiel de centre de masse. 
\end{definition}
Notons $\boldsymbol r_1=\boldsymbol{GA_1}$ et $\boldsymbol r_2=\boldsymbol{GA_2}$ les rayons-vecteurs des astres dans $\mathcal R$~:
\begin{align*}
&m_1\boldsymbol r_1+m_2\boldsymbol r_2=\boldsymbol 0\\
&\boldsymbol r_2=-\frac{m_1}{m_2}\,\boldsymbol r_1
\end{align*}
Si l'on désigne par $S$ le Soleil et par $P$ la planète,
le centre d'inertie se trouve sur le segment de droite $[SP]$ qui relie les deux astres.
\subsection{Masse réduite}
$\alpha$ étant un réel, la RFD a la propriété suivante~:
\begin{align*}
\boldsymbol F^{\,ext}&=m\,\frac{d^2\boldsymbol r}{dt^2}\\
&=\frac{m}{\alpha}\,\frac{d^2(\alpha\boldsymbol r)}{dt^2}
\end{align*}
On pose $M=m/\alpha$ et $\boldsymbol{\mathrm r}=\alpha\boldsymbol r$~:
\begin{align*}
\boldsymbol F^{\,ext}&=M\,\frac{d^2\boldsymbol{\mathrm r}}{dt^2}
\end{align*}
où $\boldsymbol r$ et $\boldsymbol{\mathrm r}$ ont même origine.

\medskip
Plaçons nous dans le ré\-fé\-ren\-tiel galiléen ayant pour centre le centre d'inertie $G$ du système binaire.
Déterminons les rayons vecteur $\boldsymbol{GA_1}$ et $\boldsymbol{GA_2}$ des astres en fonction des données du problème,
c'est à dire en fonction du vecteur $\boldsymbol{A_2A_1}$ reliant les deux astres.
Nous avons~:
\begin{align*}
\boldsymbol{A_2A_1}&=\boldsymbol{A_2G}+\boldsymbol{GA_1}\\
&=\boldsymbol{GA_1}-\boldsymbol{GA_2}
\end{align*}
En notant $\boldsymbol r=\boldsymbol{A_2A_1}$~:
\begin{equation*}
\begin{aligned}
\boldsymbol r&=\boldsymbol r_1-\boldsymbol r_2\\
&=\boldsymbol r_1+\frac{m_1}{m_2}\,\boldsymbol r_1\\
&=\frac{m_1+m_2}{m_2}\,\boldsymbol r_1
\end{aligned}
\qquad et \qquad
\begin{aligned}
\boldsymbol r&=\boldsymbol r_1-\boldsymbol r_2\\
&=-\frac{m_2}{m_1}\,\boldsymbol r_2-\boldsymbol r_2\\
&=-\frac{m_1+m_2}{m_1}\,\boldsymbol r_1
\end{aligned}
\end{equation*}
et l'on a~:
\begin{subequations}\label{rho1et2}
\begin{empheq}[box=\greybox]{align}
&\boldsymbol r_1=\frac{m_2}{m_1+m_2}\,\boldsymbol r\\
&\boldsymbol r_2=\frac{-m_1}{m_1+m_2}\,\boldsymbol r
\end{empheq}
\end{subequations}
Le système formé par les deux astres en interaction gravitationnelle est supposé isolé,
on néglige les interactions avec l'extérieur.
\'Ecrivons la relation fondamentale de la dynamique dans le ré\-fé\-ren\-tiel galiléen de centre $G$,
où la somme des forces se réduit à la seule force d'intéraction gravitationnelle~:
\begin{equation*}
\begin{aligned}
\boldsymbol F^{\,ext}_1&=m_1\,\frac{d^2\boldsymbol r_1}{dt^2}\\
&=\frac{m_1m_2}{m_1+m_2}\frac{d^2\boldsymbol r}{dt^2}
\end{aligned}
\qquad et \qquad
\begin{aligned}
\boldsymbol F^{\,ext}_2&=m_2\,\frac{d^2\boldsymbol r_2}{dt^2}\\
&=\frac{-m_1m_2}{m_1+m_2}\frac{d^2\boldsymbol r}{dt^2}
\end{aligned}
\end{equation*}
Le signe plus ou moins dépend de la convention choisie pour définir $\boldsymbol r$,
c'est à dire $\boldsymbol r=\boldsymbol r_1-\boldsymbol r_2$ ou $\boldsymbol r=\boldsymbol r_2-\boldsymbol r_1$.
\begin{definition}
La masse réduite est telle que
\begin{empheq}[box=\greybox]{align*}
\mu&\triangleq\frac{m_1m_2}{m_1+m_2}
\end{empheq}
\end{definition}
Avec cette définition~:
\begin{align*}
\boldsymbol F^{\,ext}_1=\mu\,\frac{d^2\boldsymbol r}{dt^2}\qquad et \qquad\boldsymbol F^{\,ext}_2=-\mu\,\frac{d^2\boldsymbol r}{dt^2}
\end{align*}
La loi de l'action-réaction donne $\boldsymbol F^{\,ext}_2=-\boldsymbol F^{\,ext}_1$ si bien que les deux astres sont gouvernés par la même
équation~:
\begin{empheq}[box=\greybox]{align*}
\boldsymbol F^{\,ext}_1=\mu\,\frac{d^2\boldsymbol r}{dt^2}
\end{empheq}
Cette équation coïncide formellement avec celle d'un point matériel de masse $\mu$
en mouvement sous l'action de la force extérieure $\boldsymbol F^{\,ext}_1$,
dont le vecteur position a pour module $\|\boldsymbol r\|=\|\boldsymbol{A_2A_1}\|$
et pour origine le centre $G$.
Nous pouvons donc passer de la variable $\boldsymbol r_i$ à la variable $\boldsymbol r$ à condition d'utiliser la masse réduite.
Dans ce cas, les trajectoires $\boldsymbol r_i(t)$ s'obtiennent ensuite avec les formules \eqref{rho1et2}.
\subsection{Loi des aires}
\subsubsection*{Première méthode}
Par symétrie du problème et isotropie de l'espace,
chaque force attractive exer\-cée par un astre sur l'autre est dirigée suivant la droite passant par les deux astres,
vers le centre d'inertie du système.
Si l'on se place dans le référentiel de centre de masse, le centre d'inertie est un point fixe.
Chaque force étant constam\-ment dirigée vers ce centre fixe, le problème est dit à \textit{force centrale}.
Dans ce référentiel, les forces n'ont de composante que selon le vecteur~$\boldsymbol e_\rho$,
et le mouvement est donc plan.
En utilisant l'accélération en coordonnées polaires \eqref{acc_pol} p.~\pageref{acc_pol},
\begin{align*}
\boldsymbol F^{\,ext}_1&=\mu\,\frac{d^2\boldsymbol r}{dt^2}\\
F_1\boldsymbol e_{\rho}&=\mu\left(\ddot{\rho}-\rho\dot\theta^2\right)\boldsymbol e_\rho
+\mu\left(2\dot\rho\dot\theta+\rho\ddot\theta\right)\boldsymbol e_\theta 
\end{align*}
qui donne le système d'équations~:
\begin{subequations}\label{fmr}
\begin{empheq}[box=\greybox]{align}
&\mu\left(\ddot{\rho}-\rho\dot\theta^2\right)=F_1\label{fmr1}\\
&2\dot\rho\dot\theta+\rho\ddot\theta=0\label{fmr2}
\end{empheq}
\end{subequations}
Considèrons la seconde relation et multiplions par $\rho$ de chaque côté de l'égalité~:
\begin{align*}
2\rho\dot\rho\dot\theta+\rho^2\ddot\theta&=0\\
\frac{d}{dt}\left(\rho^2\dot\theta\right)&=0
\end{align*}
\begin{theoreme}Loi des aires\\
La relation
\begin{empheq}[box=\greybox]{align}\label{C}
\rho^2\dot\theta&=C
\end{empheq}
est appellée loi des aires, $C$ est appelée constante des aires.
\end{theoreme}
Elle s'applique à chacun des astres du système binaire.
Elle ne dépend pas de la force, elle est donc valable pour toute force centrale et pas seulement pour la force gravitationnelle.
Nous verrons que selon la \textit{première loi de Kepler},
l'orbite des planètes est une ellipse dont l'un des foyer est le centre de gravité $G$ du système planète-Soleil.
La constante des aires est alors simplement le double de l'aire balayée par unité de temps par le rayon de la planète,
ce qui justifie son nom.

\begin{figure}[H]
\begin{pspicture}(-4,-2.5)(4,3)%aire reservee a l'image
%\psgrid(-4,-2.5)(4,3)
\psline{->}(0,0)(4,0)
\psline{->}(0,0)(0,3)
\rput(4.3,0){$x$}
\rput(0,3.3){$y$}
\rput(-.2,-.2){$o$}
\rput(2.236,.3){$G$}
\rput(.5,-.2){$\boldsymbol i$}
\rput(-.2,.5){$\boldsymbol j$}
\psline{->}(0,0)(1,0)
\psline{->}(0,0)(0,1)
\SpecialCoor
\psellipse(0,0)(3,2)
\psline(2.236;0)(2;-80)
\psline(2.236;0)(2.42;-40)
\psline(2.236;0)(2.42;-140)
\psline(2.236;0)(2.15;-120)
\rput(1;-120){$\rho$}
\rput(2.8;-130){$\rho d\theta$}
\rput(1.8;-123){$s_{1}$}
\rput(2;-40){$s_{2}$}
\end{pspicture}
\caption{Aires balayées en des temps égaux}
\label{aire}
\end{figure}

Soit $d\mathcal{A}$ l'aire infinitésimale balayée par le rayon $\rho$ d'un astre en un intervalle de temps infinitésimal $dt$.
Cette aire est égale au demi-rectangle de côtés $\rho$ et $\rho d\theta$ :
\begin{align*}
d\mathcal{A}&=\tfrac{1}{2}\,\rho\times\rho d\theta\\
\frac{d\mathcal{A}}{dt}&=\tfrac{1}{2}\,\rho^2\dot\theta=\tfrac{1}{2}\,C\\
\int_{0}^{s}d\mathcal{A}&=\tfrac{1}{2}\,C\int_{0}^{\tau} dt
\end{align*}
\begin{empheq}[box=\greybox]{align*}
s&=\tfrac{1}{2}\,C\tau
\end{empheq}
Les aires $s_{1}$ et $s_{2}$ de la figure~\ref{aire} balayées en des temps égaux sont égales.
On peut énoncer la \textit{deuxième loi de Kepler} :
le rayon d'un astre balaye des surfaces d'aires égales en des temps égaux.
En une période $T$, le rayon parcours la surface totale\footnote{Voir Ellipse.pdf}, d'aire $\pi\,ab$, de l'ellipse :
\begin{align*}
\int_{0}^{\pi ab}d\mathcal{A}&=\tfrac{1}{2}\,C\int_{0}^{T} dt
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{per}
C&=\frac{2\pi ab}{T}
\end{empheq}
En introduisant la pulsation $\omega$, telle que ${\omega=2\pi/T}$,
\begin{align*}
C&=\omega\,ab
\end{align*}
\subsubsection*{Deuxième méthode}
En utilisant la vitesse en coordonnées polaires \eqref{vit_pol} p.~\pageref{vit_pol},
le vecteur moment cinétique $\boldsymbol L$ s'écrit~: 
\begin{align*}
\boldsymbol L&=\boldsymbol r\times \mu\boldsymbol v\\
&=\rho\boldsymbol e_{\rho}\times \mu(\dot\rho\boldsymbol e_\rho+\rho\dot\theta\boldsymbol e_\theta )\\
&=
\begin{pmatrix}
\rho\\
0\\
0
\end{pmatrix}
\times\mu
\begin{pmatrix}
\dot\rho\\
\rho\dot\theta\\
0
\end{pmatrix}\\
&=\mu\rho^2\dot\theta\,\boldsymbol{e}_{z}
\end{align*}
Sa dérivée par rapport au temps est égale à la somme des moments des forces exercées sur cet astre~:
\begin{align*}
\frac{d\boldsymbol L}{dt}&=\rho\boldsymbol e_{\rho}\times F_1\boldsymbol e_{\rho}\\
&=\boldsymbol 0
\end{align*}
Le moment cinétique est donc constant~:
\begin{align*}
\mu\rho^2\dot\theta\,\boldsymbol{e}_{z}&=\boldsymbol C\\
\rho^2\dot\theta&=C
\end{align*}
La trajectoire de chaque astre est donc une courbe plane, dont le plan est perpendiculaire au vecteur moment cinétique $\boldsymbol L_i$ constant
et contient les vecteurs $\boldsymbol e_{\rho}$ et $\boldsymbol e_{\theta}$.

\medskip
Servons-nous de la constante des aires ou de la norme du moment cinétique pour réécrire la relation~\eqref{fmr1} p.~\pageref{fmr1}~:
\begin{align*}
\mu\ddot{\rho}&=\mu\rho\dot\theta^2+F_1\\
&=\frac{L^2}{\mu\rho^3}+F_1
\end{align*}
Si la force $F_1$ est conservative, c'est à dire si elle dérive d'une énergie potentielle $E_p(\rho_1)$, alors~:
\begin{align*}
\mu\ddot{\rho}&=\frac{L^2}{\mu\rho^3}-\frac{dE_p(\rho_1)}{d\rho_1}\\
m_1\ddot{\rho}_1&=\frac{L^2_1}{m_1\rho^3_1}-\frac{dE_p(\rho_1)}{d\rho_1}\\
&=-\frac{d}{d\rho_1}\left(\frac{L^2_1}{2m_1\rho^2_1}+E_p(\rho_1)\right)
\end{align*}
Le mouvement radial est similaire à un mouvement à une dimension (avec $\rho_1>0$) dans un potentiel $L^2_1/(2m_1\rho^2_1)+E_p(\rho_1)$.
\begin{definition}
On appelle potentiel effectif, le potentiel~:
\begin{empheq}[box=\greybox]{align*}
V_{eff}(\rho)&=\frac{L^2}{2m\rho^2}+E_p(\rho)
\end{empheq}
\end{definition}
\section{Force de gravitation}\label{fg}
Reprenons la relation~\eqref{fmr1} p.~\pageref{fmr1} valable pour toute force centrale~:
\begin{align*}
&F_1=\mu\left(\ddot{\rho}-\rho\dot\theta^2\right)
\end{align*}
Le modèle de force pour la force de gravitation est en raison inverse du carré de la distance des deux astres.
Soit $\mathcal G$ la constante de gravitation,
\begin{align*}
-\mathcal G\,\frac{m_1m_2}{\rho^2}&=\frac{m_1m_2}{m_1+m_2}\left(\ddot{\rho}-\rho\dot\theta^2\right)
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{A}
\rho^2\ddot{\rho}-\rho^{3}\dot\theta^2&=-\mathcal G(m_1+m_2)
\end{empheq}
Cette équation différentielle du second ordre par rapport au temps est difficile à intégrer.
En raisonnant sur l'énergie plutôt que sur les forces,
nous obtiendrons la même équation différentielle intégrée une première fois par rapport au temps.
En effet, la relation fondamentale de la dynamique est équivalente à la conservation de l'énergie mécanique
lorsque les forces dérivent d'un potentiel.
Le théorème de l'énergie cinétique permet de travailler sur la vitesse
par l'intermédiaire de l'énergie cinétique, plutôt que sur l'accélération.
\section{Conservation de l'énergie mécanique}
Cherchons le potentiel de la force gravitationnelle exercée sur un astre,
autrement dit l'é\-ner\-gie po\-ten\-tiel\-le de pesanteur $E_{pp}(\rho)$ de cet astre~:
\begin{align*}
-\mathcal G\,\frac{m_1m_2}{\rho^2}\,\boldsymbol e_\rho&=-\grad E_{pp}(\rho)\\
\mathcal G\,\frac{m_1m_2}{\rho^2}\,\boldsymbol e_\rho&=\frac{\partial E_{pp}(\rho)}{\partial\rho}\,\boldsymbol e_\rho\\
\mathcal G\,\frac{m_1m_2}{\rho^2}&=\frac{dE_{pp}(\rho)}{d\rho}
\end{align*}
Seule une différence d'énergie potentielle intervient si bien qu'elle est toujours définie à une constante additive près,
que l'on ajoute puis que l'on soustrait.
Ajuster la valeur de cette constante nous permet de définir le zéro de l'énergie potentielle où cela nous convient le mieux.
En annulant l'énergie potentielle à l'infini~:
\begin{align*}
dE_{pp}&=\mathcal G\,\frac{m_1m_2}{\rho^2}\,d\rho\\
\int_{\rho=+\infty}^{\rho}dE_{pp}&=\mathcal Gm_1m_2\int_{\rho=+\infty}^{\rho}\frac{d\rho}{\rho^2}\\
\left[E_{pp}\right]_{\rho=+\infty}^{\rho}&=-\mathcal Gm_1m_2\left[\frac{1}{\rho}\right]_{\rho=+\infty}^{\rho}
\end{align*}
\begin{empheq}[box=\greybox]{align*}
E_{pp}&=-\mathcal G\,\frac{m_1m_2}{\rho}
\end{empheq}
Nous pouvons aussi raisonner sur la constante d'intégration~:
\begin{align*}
\int dE_{pp}&=\mathcal Gm_1m_2\int\frac{d\rho}{\rho^2}\\
E_{pp}&=-\mathcal G\,\frac{m_1m_2}{\rho}+C^{ste}
\end{align*}
Prendre une constante nulle revient à poser $E_{pp}=0$ en ${\rho=+\infty}$.

\medskip
\'Ecrivons maintenant l'expression de l'énergie cinétique~:
\begin{align*}
E_c&=\tfrac{1}{2}\mu\boldsymbol v^2
\end{align*}
\begin{empheq}[box=\greybox]{align*}
E_c&=\tfrac{1}{2}\mu\left(\dot{\rho}^2+\rho^2\dot\theta^2\right)
\end{empheq}
\begin{rmq}
Nous pouvons aussi prendre le potentiel effectif,
\begin{align*}
V_{eff}(\rho)&=\frac{L^2(\mu)}{2\mu\rho^2}+E_p(\rho)\\
&=\tfrac{1}{2}\mu\rho^2\dot\theta^2+E_p(\rho)
\end{align*}
et l'énergie cinétique~:
\begin{align*}
E_c&=\tfrac{1}{2}\mu\dot{\rho}^2 
\end{align*}
\end{rmq}
La conservation de l'énergie mécanique donne l'équation différentielle du mouvement intégrée une première fois~:
\begin{align*}
&E_c+E_{pp}=E\\
&\frac{m_1m_2}{2(m_1+m_2)}\left(\dot{\rho}^2+\rho^2\dot\theta^2\right)-\mathcal G\,\frac{m_1m_2}{\rho}=E
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{equadiff}
\dot{\rho}^2+\rho^2\dot\theta^2-\frac{2\mathcal G(m_1+m_2)}{\rho}&=\frac{2E}{\mu}
\end{empheq}
qui est du premier ordre par rapport au temps.

\medskip
En la dérivant par rapport au temps on vérifie qu'elle redonne l'équation différentielle~\eqref{A}~:
\begin{align*}
\frac{d}{dt}\left[\dot{\rho}^2+\rho^2\dot\theta^2-\frac{2\mathcal G(m_1+m_2)}{\rho}\right] &=\frac{d}{dt}\left(\frac{2E}{\mu}\right) \\
2\dot\rho\ddot{\rho}+2\rho\dot\rho\dot\theta^2+2\rho^2\dot\theta\ddot\theta+2\mathcal G(m_1+m_2)\frac{\dot\rho}{\rho^2}&=0
\end{align*}
En dérivant la loi des aires~:
\begin{align*}
\frac{d}{dt}\,\left(\rho^2\dot\theta\right)&=\frac{d}{dt}\,C\\
2\rho\dot\rho\dot\theta+\rho^2\ddot\theta&=0\\
\rho^2\ddot\theta\times2\dot\theta&=-2\rho\dot\rho\dot\theta\times2\dot\theta\\
2\rho^2\dot\theta\ddot\theta&=-4\rho\dot\rho\dot\theta^2
\end{align*}
Par conséquent~:
\begin{align*}
2\dot\rho\ddot{\rho}-2\rho\dot\rho\dot\theta^2+2\mathcal G(m_1+m_2)\frac{\dot\rho}{\rho^2}=0&\\
\ddot{\rho}-\rho\dot\theta^2+\mathcal G(m_1+m_2)\frac{1}{\rho^2}=0&\\
\rho^2\ddot{\rho}-\rho^{3}\dot\theta^2=-\mathcal G(m_1+m_2)&
\end{align*}
\section{Résolution de l'équation différentielle du mouvement}
Revenons à l'équation différentielle \eqref{equadiff} p.~\pageref{equadiff} du premier ordre par rapport au temps~:
\begin{align*}
\dot{\rho}^2+\rho^2\dot\theta^2-\frac{2\mathcal G(m_1+m_2)}{\rho}&=\frac{2E}{\mu}
\end{align*}
Si nous voulons l'équation de la trajectoire de la forme ${\rho=f(\theta)}$,
il faut supprimer le temps dans cette équation différentielle.
Pour cela utilisons la constante des aires \eqref{C} p.~\pageref{C} : 
\begin{align*}
\rho^2\frac{d\theta}{dt}&=C
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{lda}
dt&=\frac{\rho^2d\theta}{C}
\end{empheq}
Le premier terme $\dot{\rho}^2$ de l'équation différentielle se transforme comme ceci~:
\begin{align*}
\dot{\rho}^2&=\left(\frac{d\rho}{dt}\right)^2\\
&=C^2\left(\frac{d\rho}{\rho^2d\theta}\right)^2
\end{align*}
Or,
\begin{align*}
\frac{d}{d\theta}\left(\rho^{-1}\right)&=-\rho^{-2}\,\frac{d\rho}{d\theta}\\
\frac{d(1/\rho)}{d\theta}&=-\frac{d\rho}{\rho^2d\theta}
\end{align*}
donc,
\begin{align*}
\dot{\rho}^2&=C^2\left[\frac{d(1/\rho)}{d\theta}\right]^2
\end{align*}
Les deux termes suivants se transforment comme ceci~:
\begin{align*}
\rho^2\dot\theta^2-\frac{2\mathcal G(m_1+m_2)}{\rho}&=\frac{C^2}{\rho^2}-\frac{2\mathcal G(m_1+m_2)}{\rho}\\
&=\left[\frac{C}{\rho}-\frac{\mathcal G(m_1+m_2)}{C}\right]^2-\left[\frac{\mathcal G(m_1+m_2)}{C}\right]^2
\end{align*}
L'équation différentielle \eqref{equadiff} p.~\pageref{equadiff} s'écrit maintenant~:
\begin{align*}
C^2\left[\frac{d(1/\rho)}{d\theta}\right]^2+\left[\frac{C}{\rho}-\frac{\mathcal G(m_1+m_2)}{C}\right]^2
&=\left[\frac{\mathcal G(m_1+m_2)}{C}\right]^2+\frac{2E}{\mu}\\
\left[\frac{d(1/\rho)}{d\theta}\right]^2+\left[\frac{1}{\rho}-\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2
&=\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}
\end{align*}
Nous n'avons plus que des termes en $1/\rho$.
Effectuons un premier changement de variable pour travailler avec l'inverse de la variable~$\rho$ :
\begin{align*}
\mathrm{r}&=\frac{1}{\rho}
\end{align*}
L'équation différentielle devient~:
\begin{align*}
\left[\frac{d\mathrm{r}}{d\theta}\right]^2+\left[\mathrm{r}-\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2
&=\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}
\end{align*}
Effectuons un deuxième changement de variable pour s'affranchir de la constante additive~:
\begin{align*}
\mathsf{r}&=\mathrm{r}-\frac{\mathcal G(m_1+m_2)}{C^2}
\end{align*}
L'équation différentielle devient~:
\begin{align*}
\left(\frac{d\mathsf{r}}{d\theta}\right)^2+\mathsf{r}^2&=\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}
\end{align*}
Effectuons un troisième et dernier changement de variable~:
\begin{align*}
r&=\frac{\mathsf{r}}{\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}}
\end{align*}
d'où,
\begin{align*}
\mathsf{r}&=r\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}\\
\mathsf{r}^2&=r^2\left\{ \left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}\right\} 
\end{align*}
et,
\begin{align*}
\frac{d\mathsf{r}}{d\theta}&=\frac{dr}{d\theta}\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}\\
\left(\frac{d\mathsf{r}}{d\theta}\right)^2
&=\left(\frac{dr}{d\theta}\right) ^2\left\{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}\right\}
\end{align*}
L'équation différentielle devient,
\begin{align*}
\left[\left(\frac{dr}{d\theta}\right) ^2+r^2\right]\left\{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}\right\}
&=\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}
\end{align*}
soit~:
\begin{empheq}[box=\greybox]{align}\label{arcos}
\left(\frac{dr}{d\theta}\right)^2+r^2&=1
\end{empheq}

\medskip
\textit{Première méthode de résolution}
\begin{align*}
\frac{dr}{d\theta}&=\pm(1-r^2)^{1/2}\\
d\theta&=\pm\frac{dr}{(1-r^2)^{1/2}}
\end{align*}
Lorsque le signe est positif cette équation s'intègre en arcosinus plus un angle constant $\phi$~:
\begin{align*}
\theta&=\arccos(r)+\phi\\
r&=\cos(\theta-\phi)
\end{align*}
Lorsque le signe est négatif cette équation s'intègre en arcsinus plus un angle constant $\beta$~:
\begin{align*}
\theta&=\arcsin(r)+\beta\\
r&=\sin(\theta-\beta)
\end{align*}
Pour $\beta=\phi-\pi/2$,
\begin{align*}
r&=\sin(\theta-\phi+\pi/2)\\
&=\cos(\theta-\phi)
\end{align*}
Les solutions étant identiques, on conserve ${r=\cos(\theta-\phi)}$.

\medskip
\textit{Deuxième méthode de résolution}\\
Dérivons par rapport à $\theta$ :
\begin{align*}
\frac{d}{d\theta}\left[\left(\frac{dr}{d\theta}\right) ^2+r^2\right]=0&\\
2\,\frac{dr}{d\theta}\,\frac{d^2r}{d\theta^2}+2\,r\,\frac{dr}{d\theta}=0&\\
\frac{d^2r}{d\theta^2}+r=0&\\
r=\cos(\theta-\phi)&
\end{align*}
Revenons à la variable $\rho$ :
\begin{align*}
\cos(\theta-\phi)&=r\\
&=\frac{\mathsf{r}}{\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}}\\
&=\frac{\mathrm{r}-\frac{\mathcal G(m_1+m_2)}{C^2}}{\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}}\\
&=\frac{\frac{1}{\rho}-\frac{\mathcal G(m_1+m_2)}{C^2}}{\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}}
\end{align*}
d'où~:
\begin{align*}
\frac{1}{\rho}&=\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}\ \cos(\theta-\phi)+\frac{\mathcal G(m_1+m_2)}{C^2}
\end{align*}
En posant $p$ le \textit{paramètre} de l'ellipse, tel que
\begin{align}
p=\frac{C^2}{\mathcal G(m_1+m_2)}\label{para}
\end{align}
et $e$ \textit{l'excentricité} de l'ellipse, telle que
\begin{align}
e=\frac{C^2}{\mathcal G(m_1+m_2)}\,\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}\label{excentricite}
\end{align}
l'équation s'écrit~:
\begin{align*}
\frac{1}{\rho}&=\frac{e}{p}\,\cos(\theta-\phi)+\frac{1}{p}
\end{align*}
\begin{empheq}[box=\greybox]{align}
\rho&=\frac{p}{1+e\cos(\theta-\phi)}\label{conique}
\end{empheq}
C'est l'équation polaire d'une conique, intersection d'un cône de ré\-vo\-lu\-tion avec un plan.
L'équation de la trajectoire du premier astre s'obtient grâce à la relation~\eqref{rho1et2} p.~\pageref{rho1et2}~:
\begin{align*}
\rho_1&=\frac{m_2}{m_1+m_2}\,\rho\\
&=\frac{m_2p}{(m_1+m_2)[1+e\cos(\theta-\phi)]} 
\end{align*}
De même, l'équation de la trajectoire du second astre s'écrit~:
\begin{align*}
\rho_2&=\frac{m_1}{m_1+m_2}\,\rho\\
&=\frac{m_1p}{(m_1+m_2)[1+e\cos(\theta-\phi)]} 
\end{align*}
\section{La trajectoire}
Dans ce qui suit nous considérons le système formé d'une planète quelconque du système solaire et du Soleil, en le supposant isolé.
\subsection{L'excentricité e}
L'excentricité \eqref{excentricite} s'écrit~:
\begin{align*}
e&=\frac{C^2}{\mathcal G(m_1+m_2)}\,\sqrt{\left[\frac{\mathcal G(m_1+m_2)}{C^2}\right]^2+\frac{2E}{\mu C^2}}
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{excen}
e&=\sqrt{1+\frac{2EC^2}{\mathcal G^2m_1m_2(m_1+m_2)}}
\end{empheq}
Elle définit la forme de l'orbite.
\begin{itemize}
\item Si $E=-\mathcal G^2m_1m_2(m_1+m_2)/\left(2C^2\right)$ alors $e=0$ et $\rho=p$~:
l'orbite est un cercle. Nous voyons que $p$ est le paramètre de taille de l'orbite
\item Si $-\mathcal G^2m_1m_2(m_1+m_2)/\left(2C^2\right)<E<0$ alors $0<e<1$ : l'orbite est une ellipse
\item Si $E=0$ alors $e=1$~: l'orbite est une parabole
\item Si $E>0$ alors $e>1$~: l'orbite est une hyperbole
\end{itemize}
Remarquons que l'excentricité est positive ou nulle.
\subsection{L'ellipse}
Pour $0<e<1$ et pour $\phi=0$, l'équation \eqref{conique} p.~\pageref{conique}
représente dans le repère $(F,x,y)$ l'équation polaire d'une ellipse centrée
sur $o(-c,0)$, dont le rayon $\rho$ est mesuré depuis l'un des foyers.

\begin{figure}[H]
\begin{pspicture}(-6,-3)(6,3.5)%aire reservee a l'image
%\psgrid(-6,-3)(6,3.5)
\psline{->}(-4.5,0)(4.5,0)
\psline{->}(2.236,-2.5)(2.236,3)
\rput(4.8,0){$x$}
\rput(2.236,3.3){$y$}
\rput(0,.3){$o$}
\rput(1.9,-.3){$F$}
\rput(-2.136,.4){$F'$}
\SpecialCoor
\psellipse(0,0)(3,2)
\psline(2.236,0)(1.07,1.853)
\rput(1.2,2.3){$M(x,y)$}
\rput(1.4,.8){$\rho$}
\psarc{->}(2.236,0){.5}{0}{123}
\rput(2.6,.6){$\theta$}
\rput(-2.236,0){+}
\psline{<->}(0,0)(2.236,0)
\rput(1.1,-.2){$c$}
\rput(-.3,-1){$b$}
\rput(-1.5,-.3){$a$}
\psline{<->}(0,0)(0,-2)
\psline{<->}(0,0)(-3,0)
\end{pspicture}
\caption{Orbite elliptique}
\label{ell}
\end{figure}

L'angle $\theta$ est appelé \textit{anomalie vraie}.
Nous pouvons énoncer la \textit{première loi de Kepler}~: les planètes dé\-cri\-vent des ellipses dont le centre de gravité
du système planète-Soleil occupe l'un des foyers.

\begin{rmq}
Le Soleil décrit également une ellipse dont le centre de gravité du système planète-Soleil occupe l'un des foyers.
\`A l'exception de Jupiter le centre de gravité du système planète-Soleil est dans le Soleil.
\end{rmq}

Le couple $(p,e)$ définit une ellipse de manière unique, sa taille et sa forme.
En coordonnées rectangulaires, l'équation d'une ellipse dont le grand axe est sur l'axe des abscisses, s'écrit~:
\begin{align*}
\frac{x^2}{a^2}+\frac{y^2}{b^2}=1
\end{align*}
Le couple formé du demi-grand axe $a$ et du demi-petit axe $b$ permet donc aussi de définir parfaitement une ellipse.
En étudiant le passage des coordonnées polaires $(\rho,\theta)$ vers les coordonnées rectangulaires $(x,y)$
de l'équation d'une même ellipse, nous obtenons les relations suivantes~:
\begin{subequations}
\begin{empheq}[box=\greybox]{align}
&p=\frac{b^2}{a}\label{pab}\\
&e=\sqrt{1-\frac{b^2}{a^2}}\label{eab}
\end{empheq}
\end{subequations}
\subsection{Orientation de l'ellipse}
\`A partir de l'équation polaire d'une conique \eqref{conique} p.~\pageref{conique}, pour $\theta=\phi+2k\pi,\ k\in\mathbb{Z}$
nous avons~:
\begin{align*}
\theta-\phi &=2k\pi\\
\cos(\theta-\phi)&=1\\
\rho&=\frac{p}{1+e}\\
&=\rho_{min}
\end{align*}
Nous sommes au périhélie, point de l'orbite le plus proche du Soleil.
Le grand axe de l'ellipse fait donc un angle $\phi$ avec l'axe des abscisses.

\begin{figure}[H]
\begin{pspicture}(-6,-3)(6,4)%aire reservee a l'image
%\psgrid(-6,-3)(6,4)
\rput{0}(2.160,.579){
\psline{->}(-1,0)(3,0)
\psline{->}(0,-2)(0,2)
\rput(3.3,0){$x$}
\rput(0,2.3){$y$}
\psarc{->}(0,0){1.5}{0}{15}
\rput(1.8;7){$\phi$}
}
\put(1.8,1.2){$\theta$}
\rput(1.1,1.2){$\rho$}
\rput{15}(0,0){
    \psline[linestyle=dotted](-3.5,0)(4.5,0)
    \psline[linestyle=dotted](0,-2.5)(0,2.5)
    \SpecialCoor
    \psellipse(0,0)(3,2)
    \psline(2.236,0)(1.07,1.853)
    \psarc{->}(2.236,0){.5}{-15}{123}
}
\end{pspicture}
\caption{Orientation de l'ellipse}
\label{elli}
\end{figure}

Pour $\theta=\phi+(2k+1)\pi,\ k\in\mathbb{Z}$ nous avons~:
\begin{align*}
\theta-\phi &=\pi+2k\pi\\
\cos(\theta-\phi)&=-1\\
\rho&=\frac{p}{1-e}\\
&=\rho_{max}
\end{align*}
Nous sommes à l'aphélie, point de l'orbite le plus éloigné du Soleil.

\medskip
L'orbite de la Terre autour du Soleil est dans un plan appelé \textit{plan de l'écliptique},
qui contient à chaque instant la Terre et le Soleil.
La révolution terrestre donne lieu au mouvement apparent annuel du Soleil à travers les constellations du zodiac.
Le plan de l'écliptique fait un angle de $\ang{23;26;}$ avec le plan de l'équateur terrestre, et leur intersection forme une droite dans l'espace.
Chaque année vers le 21 septembre, le Soleil passe à travers le plan de l'équateur terrestre en allant du Nord au Sud.
On pose $\theta=\ang 0$ à ce moment précis appelé \textit{équinoxe d'automne} (EA), en retour cela fixe la valeur $\phi\approx\ang{103}$.
Sur la figure \ref{elli2} l'axe des abscisses est confondu avec la droite d'intersection des deux plans.
Dans l'hémisphère Nord, le Soleil de midi monte de moins en moins haut dans le ciel et en $\theta=\pi/2$ vers le 21 décembre il atteint son point le plus bas dans le ciel, au \textit{solstice d'hiver} (SH). Le Soleil de midi remonte alors dans le ciel et en $\theta=\pi$ vers le 21 mars il passe à nouveau à travers le plan de l'équateur terrestre, cette fois en allant du Sud au Nord, à \textit{l'équinoxe de printemps} (EP).
La droite passant par les équinoxes est donc la droite d'intersection des deux plans.
Le Soleil de midi continue de monter et en $\theta=3\pi/2$ vers le 21 juin il atteint son point le plus haut dans le ciel, au \textit{solstice d'été} (SE).

\begin{figure}[H]
\begin{pspicture}(-6,-4)(6,5)%aire reservee a l'image
%\psgrid(-6,-4)(6,5)
\rput{0}(-.503,2.179){
    \psline{->}(-2,0)(3,0)
    \psline{->}(0,-6)(0,2)
    \rput(3.3,0){$x$}
    \rput(0,2.3){$y$}
    \rput(-1.5,.3){$EP$}
    \rput(2,.3){$EA$}
    \rput(.4,1.2){$SH$}
    \rput(-.4,-5.2){$SE$}
    \psarc{->}(0,0){.5}{0}{103}
    \rput(.7;50){$\phi$}
}
\rput{103}(0,0){
    \psline[linestyle=dotted](-3.5,0)(4.5,0)
    \psline[linestyle=dotted](0,-2.5)(0,2.5)
    \SpecialCoor
    \psellipse(0,0)(3,2)
}
\end{pspicture}
\caption{Orientation de l'ellipse terrestre}
\label{elli2}
\end{figure}

\subsection{Demi-grand axe a}\label{a}
Soit $2a$ la longueur du grand axe de l'ellipse~:
\begin{align*}
a&=\frac{\rho_{min}+\rho_{max}}{2}\\
&=\frac{1}{2}\left(\frac{p}{1+e}+\frac{p}{1-e}\right)\\
&=\frac{p}{1-e^2}
\end{align*}
en utilisant les relations~\eqref{para} p.~\pageref{para} et \eqref{excen} p.~\pageref{excen}~:
\begin{align}
a&=\frac{C^2}{\mathcal G(m_1+m_2)}\times\left(-\frac{\mathcal G^2m_1m_2(m_1+m_2)}{2EC^2}\right)\notag\\
&=-\frac{\mathcal Gm_1m_2}{2E}\label{grand}
\end{align}
On note au passage la relation
\begin{align}\label{pa}
p&=a(1-e^2)
\end{align}
Nous pouvons exprimer le rayon minimum en fonction du demi-grand axe $a$,
\begin{align*}
\rho_{min}&=\frac{p}{1+e}\\
&=\frac{a(1-e^2)}{1+e}\\
&=a(1-e)\\
&=a-ae
\end{align*}
de même pour le rayon maximum
\begin{align*}
\rho_{max}&=\frac{p}{1-e}\\
&=\frac{a(1-e^2)}{1-e}\\
&=a(1+e)\\
&=a+ae
\end{align*}
\subsection{Demi-petit axe b}\label{b}
En reprenant la relation \eqref{eab} p.~\pageref{eab},
\begin{align*}
e&=\sqrt{1-\frac{b^2}{a^2}}\\
e^2-1&=-\frac{b^2}{a^2}\\
b&=a\sqrt{1-e^2}
\end{align*}
en utilisant les relations \eqref{excen} p.~\pageref{excen} et \eqref{grand} p.~\pageref{grand},
\begin{align}
b&=a\sqrt{\frac{-2EC^2}{\mathcal G^2m_1m_2(m_1+m_2)}}\notag\\
&=aC\sqrt{\frac{1}{a\mathcal G(m_1+m_2)}}\notag\\
&=C\sqrt{\frac{a}{\mathcal G(m_1+m_2)}}\label{dpab}
\end{align}
\subsection{Période de révolution T}\label{periode}
A partir des équations \eqref{per} p.~\pageref{per} et \eqref{dpab} p.~\pageref{dpab},
et des paragraphes \ref{a} et \ref{b} p.~\pageref{b}, nous pouvons évaluer la période de révolution en fonction du demi-grand axe $a$ seul~:
\begin{align*}
T&=\frac{2\pi ab}{C}
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{T}
T^2&=\frac{4\pi^2\,a^{3}}{\mathcal G(m_1+m_2)}
\end{empheq}
Ceci constitue la \textit{troisième loi de Kepler}~:
le carré de la période de révolution est proportionnel au cube du demi-grand axe de l'ellipse.
\section{Equations paramétriques}
\subsection{L'angle en fonction du temps}
Les équations paramétriques du temps vont nous permettre de connaître la position de la planète en fonction du temps.\\
Cherchons tout d'abord ${\theta=\theta(t)}$.
En utilisant la loi des aires~\eqref{lda} p.~\pageref{lda} et l'équation de l'orbite \eqref{conique} p.~\pageref{conique} :
\begin{align*}
dt&=\frac{\rho^2}{C}\,d\theta\\
dt&=\frac{p^2d\theta}{C[1+e\cos(\theta-\phi)]^2}\\
\int_{0}^{\tau}dt&=\frac{p^2}{C}\int_{\phi}^{\vartheta}\frac{d\theta}{[1+e\cos(\theta-\phi)]^2}
\end{align*}
où l'on a pris comme condition initiale que la planète est au périhélie, c'est à dire $\theta=\phi$ à $t=0$.
Le membre de droite s'intègre en termes de fonctions élémentaires. Le résultat est trouvé grâce au logiciel de calcul formel Maxima.
Pour une trajectoire elliptique, $0<e<1$, nous avons~:
\begin{align*}
\tau&=\frac{2}{(1-e^2)^{3/2}}\arctan\left \{\frac{(1-e)\sin (\vartheta-\phi)}{\sqrt{1-e^2}[\cos (\vartheta-\phi)+1]}\right \}\\
&\qquad\qquad\qquad-\frac{2e\sin(\vartheta-\phi)}{\frac{(e^{3}-e^2-e+1)\sin^2(\vartheta-\phi)}{[\cos (\vartheta-\phi)+1]}
-(e^{3}+e^2-e-1)[\cos(\vartheta-\phi)+1]}
\end{align*}
Cette re\-la\-ti\-on est trop complexe pour que son in\-ver\-si\-on, ${\theta=\theta(t)}$, ne soit possible analytiquement.
\subsection{Le rayon en fonction du temps}
Cherchons à présent ${\rho=\rho(t)}$ en revenant à l'équation différentielle \eqref{equadiff} p.~\pageref{equadiff}~:
\begin{align*}
\dot{\rho}^2+\rho^2\dot\theta^2-\frac{2\mathcal G(m_1+m_2)}{\rho}&=\frac{2E}{\mu}\\
\left(\frac{d\rho}{dt}\right)^2+\left(\frac{C}{\rho}\right)^2-\frac{2\mathcal G(m_1+m_2)}{\rho}&=\frac{2E}{\mu}\\
\frac{d\rho}{dt}&=\pm\sqrt{\frac{2E}{\mu}+\frac{2\mathcal G(m_1+m_2)}{\rho}-\left(\frac{C}{\rho}\right)^2}\\
dt&=\frac{\pm\rho d\rho}{\sqrt{\frac{2E}{\mu}\rho^2+2\mathcal G(m_1+m_2)\rho-C^2}}
\end{align*}
A partir des relations~\eqref{para} p.~\pageref{para} et \eqref{pa} p.~\pageref{pa},
\begin{align*}
C^2&=\mathcal G(m_1+m_2)p\\
&=\mathcal G(m_1+m_2)a(1-e^2)
\end{align*}
et avec la relation \eqref{excen} p.~\pageref{excen}~:
\begin{align*}
1+\frac{2EC^2}{\mathcal G^2m_1m_2(m_1+m_2)}&=e^2\\
\frac{2E}{m_1m_2}&=\frac{(e^2-1)\mathcal G^2(m_1+m_2)}{C^2}\\
&=\frac{\mathcal G^2(m_1+m_2)(e^2-1)}{\mathcal G(m_1+m_2)a(1-e^2)}\\
\frac{2E}{\mu}&=-\frac{\mathcal G(m_1+m_2)}{a}
\end{align*}
D'où,
\begin{align}
dt&=\frac{\pm\rho d\rho}{\sqrt{-\frac{\mathcal G(m_1+m_2)}{a}\rho^2+2\mathcal G(m_1+m_2)\rho-\mathcal G(m_1+m_2)a(1-e^2)}}\notag\\
&=\pm\sqrt{\frac{a}{\mathcal G(m_1+m_2)}}\,\frac{\rho d\rho}{\sqrt{-\rho^2+2a\rho-a^2(1-e^2)}}\notag\\
&=\sqrt{\frac{a}{\mathcal G(m_1+m_2)}}\frac{\rho d\rho}{\sqrt{a^2e^2-(a-\rho)^2}}\label{dt}
\end{align}
où l'on garde uniquement le signe positif car $dt$ est positif.
Nous avons vu au paragraphe \ref{a} p.~\pageref{a} que le rayon $\rho$ varie entre $a-ae$ et $a+ae$.
Ceci suggère d'effectuer le changement de variable suivant~:
\begin{empheq}[box=\greybox]{align}\label{rho}
\rho=a-ae\cos \psi,\ \psi\in[0,\pi]
\end{empheq}
L'angle $\psi$, appelé \textit{anomalie excentrique} remplace le rayon vecteur.
Nous avons alors,
\begin{align*}
\frac{d\rho}{d\psi}&=ae\sin \psi\\
d\rho&=ae\sin \psi\,d\psi
\end{align*}
Avec cette nouvelle variable~:
\begin{align*}
dt&=\sqrt{\frac{a}{\mathcal G(m_1+m_2)}}\,\frac{(a-ae\cos \psi)\,ae\sin \psi\,d\psi}{\sqrt{a^2e^2-a^2e^2\cos^2\psi}}\\
&=\sqrt{\frac{a}{\mathcal G(m_1+m_2)}}\,\frac{(a-ae\cos \psi)\,ae\sin \psi\,d\psi}{ae\sin \psi}\\
&=\frac{a^{3/2}}{\sqrt{\mathcal G(m_1+m_2)}}\,(1-e\cos \psi)\,d\psi
\end{align*}
On suppose de nouveau qu'à l'instant initial la planète est au périhélie,
donc à l'instant initial ${\rho=\rho_{min}}$, $\psi=0$ et $\theta=\phi$.
\begin{align}
\int_{0}^{\tau} dt&=\frac{a^{3/2}}{\sqrt{\mathcal G(m_1+m_2)}}\int_{0}^{\psi}(1-e\cos \psi)\,d\psi\notag\\
\tau&=\frac{a^{3/2}}{\sqrt{\mathcal G(m_1+m_2)}}\Big[\psi-e\sin \psi\Big]_{0}^{\psi}\notag\\
&=\frac{a^{3/2}}{\sqrt{\mathcal G(m_1+m_2)}}\,(\psi-e\sin\psi)\label{tau}
\end{align}
La planète effectue une demi-période $T/2$ quand $\rho$ varie de $\rho_{min}$ à $\rho_{max}$,
et avec le changement de variable \eqref{rho} p.~\pageref{rho},
quand $\psi$ varie de $0$ à $\pi$~:
\begin{align*}
\frac{T}{2}&=\frac{a^{3/2}}{\sqrt{\mathcal G(m_1+m_2)}}\,(\pi-e\,sin\pi)\\
T&=\frac{2\pi a^{3/2}}\,{\sqrt{\mathcal G(m_1+m_2)}}
\end{align*}
On retrouve bien la période de révolution \eqref{T} p.~\pageref{T}.

\medskip
En introduisant la période $T$ puis la pulsation $\omega$ dans la relation \eqref{tau},
la solution de l'é\-qua\-ti\-on différentielle devient~:
\begin{align*}
\tau&=\frac{T}{2\pi}(\psi-e\,\sin\psi)
\end{align*}
\begin{empheq}[box=\greybox]{align}\label{eqkepler}
\omega\,\tau&=\psi-e\,\sin\psi
\end{empheq}
appelée \textit{équation de Kepler}. L'angle~$\omega\,\tau$ s'appelle \textit{a\-no\-ma\-lie moy\-en\-ne}.

Revenons à la variable~$\rho$ en inversant le changement de variable \eqref{rho} p.~\pageref{rho}.
Dans la pratique cela n'est jamais fait, nous allons voir pourquoi.
\begin{align*}
\rho&=a-ae\cos \psi\\
\frac{a-\rho}{ae}&=\cos\psi\\
\psi&=\arccos\left(\frac{a-\rho}{ae}\right) 
\end{align*}
ce qui donne
\begin{align}
\omega\tau&=\arccos\left(\frac{a-\rho}{ae}\right)-e\,\sin\arccos\left(\frac{a-\rho}{ae}\right)\label{sinarcos} 
\end{align}
On fixe $\rho$ pour obtenir $\tau$.
Pour la détermination de l'angle~$\theta$ à partir d'un rayon $\rho$ fixé,
nous utilisons l'équation de l'orbite \eqref{conique} p.~\pageref{conique}~:
\begin{align*}
\rho+\rho\,e\cos(\theta-\phi)&=p\\
\cos(\theta-\phi)&=\frac{p-\rho}{\rho e}\\
\theta&=\arccos\left(\frac{p-\rho}{\rho e}\right)+\phi 
\end{align*}
Cette méthode n'est pas optimale pour un calcul sur ordinateur, avec un sinus d'arcosinus dans l'expression \eqref{sinarcos}. 
Il est préférable de fixer $\psi$ pour trouver l'époque $\tau$ grâce à l'équation de Kepler~\eqref{eqkepler},
puis de trouver la distance au Soleil $\rho$ en utilisant le changement de variable \eqref{rho} p.~\pageref{rho}.
Enfin, nous pouvons trouver la valeur correspondante de l'angle $\theta$ en utilisant
l'équation de l'orbite \eqref{conique} p.~\pageref{conique} et le changement de variable \eqref{rho} qui donne chacun le rayon~$\rho$~:
\begin{align*}
a-ae\cos \psi&=\frac{p}{1+e\cos(\theta-\phi)}\\
1+e\cos(\theta-\phi)&=\frac{a(1-e^2)}{a(1-e\cos \psi)}\\
e\cos(\theta-\phi)&=\frac{1-e^2-(1-e\cos \psi)}{1-e\cos \psi}\\
\cos(\theta-\phi)&=\frac{\cos \psi-e}{1-e\cos \psi}\\
\theta&=\arccos\left(\frac{\cos \psi-e}{1-e\cos \psi}\right) +\phi
\end{align*}
Nous pouvons utiliser une forme alternative pour exprimer $\theta$ en fonction de $\psi$ :
\begin{align*}
1-\cos(\theta-\phi)&=1-\frac{\cos \psi-e}{1-e\cos \psi}\\
2\sin^2\left(\frac{\theta-\phi}{2}\right)&=\frac{1-e\cos \psi-\cos \psi+e}{1-e\cos \psi}\\
&=\frac{(1+e)(1-\cos \psi)}{1-e\cos \psi}
\end{align*}
et,
\begin{align*}
1+\cos(\theta-\phi)&=1+\frac{\cos \psi-e}{1-e\cos \psi}\\
2\cos^2\left(\frac{\theta-\phi}{2}\right)&=\frac{1-e\cos \psi+\cos \psi-e}{1-e\cos \psi}\\
&=\frac{(1-e)(1+\cos \psi)}{1-e\cos \psi}
\end{align*}
en faisant le rapport,
\begin{align*}
\tan^2\left(\frac{\theta-\phi}{2}\right)&=\frac{(1+e)(1-\cos \psi)}{(1-e)(1+\cos \psi)}\\
\tan\left(\frac{\theta-\phi}{2}\right)&=\sqrt{\frac{(1+e)}{(1-e)}}\tan\left(\frac{\psi}{2}\right)\\
\theta&=2\arctan\left[\sqrt{\frac{1+e}{1-e}}\tan\left(\frac{\psi}{2}\right)\right]+\phi
\end{align*}
\end{document}
