
\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}%
\usepackage[garamond]{mathdesign}
\usepackage{pst-eqdf,pst-node,pst-tools,pst-grad}
\usepackage{array,amsmath}
\usepackage{animate}
\usepackage{wasysym}
\usepackage{url}
\newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.1,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 0}}}
\makeatletter
%% adapté de \psRandom du package pstricks-add
%% pour rendre aléatoire la taille des étoiles
%% Manuel Luque
\newdimen\pssizeStar
\def\psset@sizeStar#1{\pssetlength\pssizeStar{#1}}
\psset@sizeStar{1pt}
\define@key[psset]{pst-eqd}{randomPoints}[1000]{\def\psk@randomPoints{#1}}
\psset[pst-eqd]{randomPoints=1000}
\define@boolkey[psset]{pst-eqd}[Pst@]{color}[true]{}
\psset[pst-eqd]{color=false}
\def\psRandomStar{\pst@object{psRandomStar}}%
\def\psRandomStar@i{\@ifnextchar({\psRandomStar@ii}{\psRandomStar@iii(0,0)(1,1)}}
\def\psRandomStar@ii(#1){\@ifnextchar({\psRandomStar@iii(#1)}{\psRandomStar@iii(0,0)(#1)}}
\def\psRandomStar@iii(#1)(#2)#3{%
  \def\pst@tempA{#3}%
  \ifx\pst@tempA\pst@empty\psclip{\psframe(#2)}\else\psclip{#3}\fi
  \pst@getcoor{#1}\pst@tempA
  \pst@getcoor{#2}\pst@tempB
  \begin@SpecialObj
  \addto@pscode{
    \pst@tempA\space /yMin exch def
    /xMin exch def
    \pst@tempB\space /yMax exch def
    /xMax exch def
    /dy yMax yMin sub def
    /dx xMax xMin sub def
    rrand srand                 % initializes the random generator
    /getRandReal { rand 2147483647 div } def
     \psk@randomPoints {
    /DS \pst@number\pssizeStar\space getRandReal mul def
    \@nameuse{psds@\psk@dotstyle}
     \ifPst@color getRandReal 1 1 sethsbcolor \fi
     getRandReal dx mul xMin add
     getRandReal dy mul yMin add
     Dot
     \ifx\psk@fillstyle\psfs@solid fill \fi stroke
    } repeat
  }%
  \end@SpecialObj
  \endpsclip
  \ignorespaces
}
\makeatother
%%%%%%%%%%%%%%%%%%
%timeline
\begin{filecontents}{kepler16Nasa.dat}
::0x0
::1
::2
::3
::4
::5
::6
::7
::8
::9
::10
::11
::12
::13
::14
::15
::16
::17
::18
::19
::20
::21
::22
::23
::24
::25
::26
::27
::28
::29
::30
::31
::32
::33
::34
::35
::36
::37
::38
::39
::40
::41
::42
::43
::44
::45
::46
::47
::48
::49
::50
::51
::52
::53
::54
::55
::56
::57
::58
::59
::60
::61
::62
::63
::64
::65
::66
::67
::68
::69
::70
::71
::72
::73
::74
::75
::76
::77
::78
::79
::80
::81
::82
::83
::84
::85
::86
::87
::88
::89
::90
::91
::92
::93
::94
::95
::96
::97
::98
::99
::100
::101
::102
::103
::104
::105
::106
::107
::108
::109
::110
::111
::112
::113
::114
::115
::116
::117
::118
::119
::120
::121
::122
::123
::124
::125
::126
::127
::128
::129
::130
::131
::132
::133
::134
::135
::136
::137
::138
::139
::140
::141
::142
::143
::144
::145
::146
::147
::148
::149
::150
::151
::152
::153
::154
::155
::156
::157
::158
::159
::160
::161
::162
::163
::164
::165
::166
::167
::168
::169
::170
::171
::172
::173
::174
::175
::176
::177
::178
::179
::180
::181
::182
::183
::184
::185
::186
::187
::188
::189
::190
::191
::192
::193
::194
::195
::196
::197
::198
::199
::200
::201
::202
::203
::204
::205
::206
::207
::208
::209
::210
::211
::212
::213
::214
::215
::216
::217
::218
::219
::220
::221
::222
::223
::224
::225
::226
::227
::228
::229
::230
::231
\end{filecontents}
%%%%%%%%%%%%%%%%%%

\title{Gravitation : une planète à deux soleils - partie 2}
\date{12 juillet 2\,012}
\begin{document}
\maketitle
\section{Explications sur la représentation à partir des données de la NASA}
Nous avons vu que le site de la NASA dédié à cette planète\footnote{\url{http://kepler.nasa.gov/Mission/discoveries/kepler16b/}}, fournit un grand nombre de renseignements sur les étoiles et la planète, ce qui permet de reconstituer leurs trajectoires respectives. Je rappelle les caractéristiques utiles pour le schéma et l'animation de ces trois corps.
\begin{center}
\begin{tabular}{|ll|}
\hline
\textbf{ Paramètres}&\textbf{Valeurs}\\ \hline
  \textit{Étoile A}&\\ \hline
  Masse, $M_A(M_{\astrosun})$&0.6897\\
  Rayon, $R_A(R_{\astrosun})$&0.6489\\ \hline
  \textit{Étoile B}&\\ \hline
  Masse, $M_B(M_{\astrosun})$&0.20255\\
  Rayon, $R_B(R_{\astrosun})$&0.22623\\ \hline
  \textit{Planète b}&\\ \hline
  Masse, $M_b(M_{\jupiter})$&0.333\\
  Rayon, $R_b(R_{\jupiter})$&0.7538\\ \hline
  \textit{Orbite de l'étoile binaire}&\\ \hline
  Période(jours)&41.076\\
  Demi-grand axe(ua) $\mathrm{a_1}$&0.22431\\
  Excentricité $\mathrm{e_1}$&0.15944\\
  Argument du périastre(deg) $\omega_1$&263.464\\ \hline
  \textit{Orbite de la planète circumbinaire} &\\ \hline
   Période(jours)&228.776\\
  Demi-grand axe(ua)$\mathrm{a_2}$&0.7048\\
  Excentricité $\mathrm{e_2}$&0.0069\\
  Argument du périastre(deg) $\omega_2$&318\\ \hline
\end{tabular}
\end{center}
Voici quelques explications qui permettront, je l'espère, de comprendre comment on peut reconstituer le mouvement de la planète ``Kepler-16b'' orbitant autour d'une étoile binaire, ainsi que ceux des 2 étoiles ``Kepler-16A'' et ``Kepler-16B'', à partir de ces données.
\subsection{Le couple d'étoiles}
Pour le couple d'étoiles $A$ et $B$, les données de la \textsc{nasa} caractérisent le mouvement dans le repère barycentrique, d'un corps fictif de masse égale à la masse réduite :
\[
\mu=\frac{m_Am_B}{m_A+m_B}
\]
placé en un point $M$ tels que $\overrightarrow{r_1}=\overrightarrow{CM}=\overrightarrow{AB}=\overrightarrow{r_B}-\overrightarrow{r_A}$, comme s'il subissait l'attraction gravitationnelle d'un corps fictif placé au barycentre, de masse égale à la somme des masses de deux étoiles\footnote{On pourra consulter le document suivant, dédié au problème des deux corps : \newline \url{http://melusine.eu.org/syracuse/pstricks/pst-eqdf/depotgit/gravitation/pb_2corps_doc.pdf}\newline
dont on trouvera aussi le fichier source à la même adresse : \newline\url{http://melusine.eu.org/syracuse/pstricks/pst-eqdf/depotgit/gravitation/}}.

Ce corps décrit une ellipse dont le centre de masse du système un foyer. Son équation paramétrique s'écrit :
\[
r_1=\frac{p_1}{1+e_1\cos(\theta-\omega_1)}
\]
On en déduit le paramètre de cette ellipse :
\[
p_1=a_1(1-e_1^2)
\]
Les équations respectives des deux étoiles se déduisent de la relation barycentrique :
\[
\overrightarrow{r_A}=-\frac{m_B}{m_A+m_B}\overrightarrow{r_1}\quad\mbox{et}\quad \overrightarrow{r_B}=\frac{m_A}{m_A+m_B}\overrightarrow{r_1}
\]
Ce sont deux ellipses de même excentricité, mais dont les paramètres et demi-grand axes sont différents :
\[
p_A=p_1\frac{m_B}{m_A+m_B}\quad ; \quad p_B=p_1\frac{m_A}{m_A+m_B}
\]
On peut en déduire les demi-grand axes :
\[
a_A=\frac{p_A}{1-e_1^2}\quad ; \quad a_B=\frac{p_B}{1-e_1^2}
\]
Les équations paramétriques s'écrivent :
\[
r_A=-\frac{p_A}{1+e_1\cos(\theta-\omega_1)}\quad ; \quad r_B=\frac{p_B}{1+e_1\cos(\theta-\omega_1)}
\]
\begin{center}
\begin{pspicture}(-5,-5)(5,5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{% pour le système binaire d'étoiles
        /a1 0.22431 def % semi-major axis
        /e1 0.15944 def % eccentricity
        /w1 263.464 def % argument of periaspe
        /MA 0.6897 def  % masse of star A
        /MB 0.20255 def % masse of star B
        /Mt MA MB add def % masse totale
        /P1 41.076 def  % period (day)
        % pour la planète
        /Mb 0.3333 def  % masse of planet b (M_Jupiter)
        /a2 0.7048 def  % semi-major axis
        /e2 0.0069 def  % eccentricity
        /w2 318 def  % argument of periaspe
        /P2 228.776 def % period (day)
% on en déduit le paramètre de chacune des ellipses
        /p1 a1 1 e1 dup mul add mul def
        /p2 a2 1 e2 dup mul add mul def
        }%
\pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
         radius w1 cos mul 5 mul
         radius w1 sin mul 5 mul){Mu}
\psdot(Mu)
\parametricplot[plotpoints=360,unit=5,linestyle=dashed]{0}{360}{
    /radius p1 1 e1 t w1 sub cos mul add div def
     radius t cos mul
     radius t sin mul}
\pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div neg def
        /radius1 radius MB Mt div mul def
         radius1 w1 cos mul 5 mul
         radius1 w1 sin mul 5 mul){MA0}
\parametricplot[plotpoints=360,unit=5,,linecolor={[rgb]{1 0.5 0}}]{0}{360}{
    /radius p1 1 e1 t w1 sub cos mul add div neg def
    /radius1 radius MB Mt div mul def
     radius1 t cos mul
     radius1 t sin mul}
\pscircle[fillstyle=solid,fillcolor=yellow!50](MA0){0.3}
\parametricplot[plotpoints=360,unit=5,linecolor=red]{0}{360}{
    /radius p1 1 e1 t w1 sub cos mul add div def
    /radius2 radius MA Mt div mul def
     radius2 t cos mul
     radius2 t sin mul}
\pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
        /radius1 radius MA Mt div mul def
        radius1 w1 cos mul 5 mul
        radius1 w1 sin mul 5 mul){MB0}
\pscircle[fillstyle=solid,fillcolor=red](MB0){0.15}
\parametricplot[plotpoints=360,unit=5,linecolor=blue]{0}{360}{
    /radius p2 1 e2 t w2 sub cos mul add div def
     radius t cos mul
     radius t sin mul}
\pnode(!/radius p2 1 e2 w2 w2 sub cos mul add div def
        radius w2 cos mul 5 mul
        radius w2 sin mul 5 mul){Mb0}
\pscircle[fillstyle=solid,fillcolor=blue](Mb0){0.07}
\end{pspicture}
\end{center}
Pour décrire complètement le mouvement des 3 corps, il faut pouvoir déterminer leurs positions respectives à un instant donné. Cela est possible puisque nous connaissons leurs périodes de révolution. Pour cela il faut résoudre \textit{l'équation de Kepler} qui est la relation entre l'\textit{anomalie moyenne} et l'\textit{anomalie excentrique}\footnote{Un angle compté à partir de la direction du périhélie porte le nom de d'\textit{anomalie} (cf. Guy Sérane)}.
%\[
%AE=AM +\mathrm{e}\sin AE
%\]

Guy Sérane dans ``\textit{Astronomie \& Ordinateur}'' publié chez Dunod en 1987 propose une intéressante démonstration géométrique de cette relation. Sur le schéma, les angles $E$  et $v$ représentent respectivement l'\textit{anomalie excentrique} et l'\textit{anomalie vraie}. Sa démonstration a été adaptée aux notations du schéma suivant :
\begin{center}
\begin{pspicture}(-8,-5)(3,5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{/par 2 def /exc 0.7 def /a_2 par 1 exc dup mul sub div def /b_2 a_2 1 exc dup mul sub sqrt mul def}%
\parametricplot[plotpoints=360,linecolor=red]{0}{360}{
    /radius par 1 exc t cos mul add div def
     radius t cos mul
     radius t sin mul}
\pnode(! a_2 exc mul neg 0){C} % centre de l'ellipse
\rput(C){\psline(! a_2 neg 0)(! a_2 0)\psline[linestyle=dashed](! 0 a_2)(! 0 a_2 neg)%
         \pscircle(0,0){!a_2}}
\pnode(0,0){O}% foyer
\pscircle*(0,0){0.05}\pscircle(0,0){0.1}
\uput[dl](C){O}\uput[d](O){F}
\pstVerb{/AV 120 def}%
\pnode(! /radius par 1 exc AV cos mul add div def
          radius AV cos mul
          radius AV sin mul){P} % planète
\uput{0.1}[ur](P){$S$}
\psdot(P)\psline(O)(P)
%\pstVerb{% équation de Kepler
%    5{/AETemp
%       AM Excentric  180 mul 3.14159 div AE sin mul add def
%      /AE AETemp def} repeat
%     }%
\pstVerb{/tan {dup sin exch cos div} def
         /arctan {
            dup 0 ge
            {1 atan}
            {neg 1 atan neg}
         ifelse
        } def
        /AE 1 exc sub 1 exc add div sqrt AV 2 div tan mul arctan 2 mul def
        /xP1 AE cos a_2 mul def
        /yP1 AE sin a_2 mul def}
\rput(C){\pnode(!xP1 yP1){P1}\pnode(!xP1 0){H}}
\rput(C){\pnode(!a_2 neg 0){A}\pnode(!0 b_2){B}\pnode(!a_2 0){Per}}
\pscustom[fillstyle=hlines]{%
    \psline(C)(Per)
    \parametricplot[plotpoints=360,linecolor=red]{0}{AV}{
    /radius par 1 exc t cos mul add div def
     radius t cos mul
     radius t sin mul}
     \psline(P)(O)}
\pscustom[fillstyle=vlines]{%
    \psline(O)(Per)
    \psarc(C){!a_2}{0}{!AE}
    \psline(P1)(F)}
\psdot(P1)\uput[u](P1){$M$}\psline(C)(P1)\psline(H)(P1)
\uput[d](H){$H$}
\psarc[linecolor=white,linewidth=0.1](O){0.7}{0}{!AV}
\psline(P)(F)
\psarc[arrowinset=0.05,arrowsize=0.1]{->}(O){0.7}{0}{!AV}
\rput(O){\uput{0.75}[!AV 2 div](0,0){\psframebox[linestyle=none,fillstyle=solid,framesep=2pt]{$v$}}}
\psarc[arrowinset=0.05,arrowsize=0.1]{->}(C){0.5}{0}{!AE}
\rput(C){\uput{0.55}[!AE 2 div](0,0){$E$}}
\pcline[offset=-0.65cm]{<->}(C)(O)
\ncput*[nrot=:U]{$c$}
\pcline{<->}(A)(C)
\ncput*[nrot=:U]{$a$}
\pcline{->}(C)(B)
\ncput*{$b$}
\rput(2.2,-0.3){périastre}
\uput[ur](Per){$P$}
\psline[arrowinset=0.05,arrowsize=0.1]{->}(2.8,-0.5)(1.5,-0.5)(Per)
\psline[arrowinset=0.05,arrowsize=0.1,linecolor=blue]{<->}(3,0)(0,0)(0,5)
\uput[r](0,4.8){$y$}\uput[u](2.9,0){$x$}
\end{pspicture}
\end{center}
$v$ est appelé \textit{anomalie vraie} et $E$ est appelé \textit{anomalie excentrique}.

On passe du cercle de centre $O$ et de rayon $a$ dit \textit{cercle auxiliaire principal}, à l'ellipse de même centre, de demi-grand axe $a$ et de demi-petit axe $b$ par une affinité de rapport $b/a$ suivant l'axe des ordonnées.
\[
\frac{HS}{HM}=\frac{b}{a}=\sqrt{1-e^2}
\]
On se place dans le repère $Fxy$, $F$ étant le foyer de l'ellipse décrite par l'astre, c'est-à dire pour le cas qui nous intéresse le barycentre des deux étoiles. Les coordonnées de l'astre $(P)$ s'expriment en fonction de $v$ ou de $E$ ; $r=FS$. On obtient facilement :
\[
\left\{
\begin{array}[m]{l}
x=r\cos v=a(\cos E-\mathrm{e})\\[1em]
y=r\sin v=a\sqrt{1-\mathrm{e}^2}\sin E\\[1em]
r=a(1-\cos E)
\end{array}
\right.
\]
Nous pouvons en déduire l'anomalie vraie $(v)$ en fonction de $E$ :
\[
\left\{
\begin{array}[m]{l}
\cos v=\dfrac{\cos E-\mathrm{e}}{1-\cos E}\\[1em]
\sin v=\dfrac{\sqrt{1-\mathrm{e}^2\sin E}}{1-\cos E}
\end{array}
\right.
\]
On utilise les formules de bissection pour calculer $v$ en fonction de $E$ :
\[
\tan^2\dfrac{v}{2}=\dfrac{1-\cos v}{1+\cos v}=\dfrac{1+\mathrm{e}}{1-\mathrm{e}}\tan^2\dfrac{E}{2}
\]
\[
\psframebox[linestyle=none,fillstyle=solid,fillcolor=cyan!20]{%
\tan\dfrac{v}{2}=\sqrt{\dfrac{1+\mathrm{e}}{1-\mathrm{e}}}\tan\dfrac{E}{2}
}
\]
\textbf{La démonstration de Guy Sérane qui établit la relation de Képler :}

Le rayon vecteur $FM$ balaye sur le cercle en partant de $P$, une aire $S_1$ telle que :

$S_1=$ (aire du secteur $POM$) $-$ (aire du triangle $OMF$) de base $OF=c=a\mathrm{e}$ et de hauteur $HM=a\sin E$ .
\[
S_1=\dfrac{E}{2}a^2-\dfrac{\mathrm{e}a^2\sin E}{2}=\frac{a^2}{2}\left(E-\mathrm{e}\sin E\right)
\]
Si $S_2$ est la surface balayée par le rayon vecteur $FS$ sur l'ellipse, pendant la même durée, en utilisant la propriété énoncée au début : ``\textit{l'ellipse se déduit du cercle principal par une affinité de rapport $b/a$ suivant l'axe des ordonnées}'', et en remarquant que la transformation ne porte que sur une dimension ($y$), on peut en déduire que les surfaces balayées sur l'ellipse et sur le cercle principal sont dans le même rapport :
\[
S_2=S_1\dfrac{b}{a}=\dfrac{ab}{2}\left(E-\mathrm{e}\sin E\right)
\]
Sur une  période de révolution du satellite(ou la planète, ou de l'astre) notée $T$, le rayon vecteur $FS$ balaye l'ellipse de surface $\pi ab$. D'après la deuxième loi de Képler, ou loi des aires, si $t-t_P$ est la durée du parcours entre le périastre et la position considérée, nous pouvons écrire :
\[
t-t_P=T\dfrac{\text{aire } S_2}{\text{aire  ellipse}}=T\dfrac{\dfrac{ab}{2}\left(E-\mathrm{e}\sin E\right)}{\pi ab}
\]
Cette méthode géométrique permet donc d'établir l'équation de Kepler :
\[
\psframebox[linestyle=none,fillstyle=solid,fillcolor=cyan!20]{%
\dfrac{2\pi}{T}(t-t_P)=E-\mathrm{e}\sin E}
\]
$\dfrac{2\pi}{T}$ est le \textit{mouvement moyen} de l'astre, et $\dfrac{2\pi}{T}(t-t_P)$ l'\textit{anomalie moyenne}.
%L'équation de Kepler :
%\[
%\psframebox[linestyle=none,fillstyle=solid,fillcolor=cyan!20]{%
%E-\mathrm{e}\sin E=\frac{2\pi}{T}(t-t_P)
%}
%\]
Cette équation permettra de connaître l'anomalie excentrique $E$ à un instant donné $(t)$. On pourra en déduire l'anomalie vraie.
\newline
Avec les données de la \textsc{Nasa}, nous pouvons maintenant déterminer les positions des 2 étoiles et de la planète dans le repère barycentrique à un instant quelconque, en résolvant l'équation de Képler.
\section{La résolution de l'équation  de Képler}
Le livre de Guy Sérane donne deux solutions : par un calcul itératif ou bien en utilisant un développement en série.
\subsection{Calcul par itérations pour obtenir l'anomalie excentrique}
On donne une valeur initiale, arbitraire, à l'anomalie excentrique, par exemple celle de l'anomalie moyenne et on effectue un nombre limité d'itérations, car la valeur se stabilise très rapidement, on peut éventuellement mettre une condition sur la précision pour sortir de la boucle. On fixe l'origine des temps à l'instant du passage au périastre $t_P=0$. Les données de la \textsc{Nasa} ne contiennent pas de renseignements sur le décalage temporel du passage au périastre des étoiles et de la planète. On suppose donc que les trois corps, à l'instant initial choisi, passent au même instant par leur périastre respectif.
\begin{verbatim}
\pstVerb{
% les angles doivent être en degrés
% pour calculer les sinus, cosinus
% en postscript
  /rad2deg { 180 mul Pi div } bind def % radians -> degrés
  /AM t 2Pi mul T div def % anomalie moyenne à la date t
  /AE AM def % valeur initiale de anomalie excentrique
% équation de Kepler
    5 {% 5 boucles
     /AETemp % anomalie excentrique provisoire
      AM excentricite rad2deg AE sin mul add
      def
     /AE AETemp def
      }repeat
       }%
\end{verbatim}
\subsection{Développement en série pour obtenir l'anomalie vraie}
Désignons par $a_m$ l'anomalie moyenne et $a_v$ l'anomalie vraie.
\[
a_v=a_m+\left(2\mathrm{e}-\dfrac{\mathrm{e}^3}{4}\right)\sin a_m +\left( \dfrac{5}{4}\mathrm{e}^2-\dfrac{11}{24}\mathrm{e}^4\right)\sin2a_m+
 \dfrac{13}{12}\mathrm{e}^3\sin3a_m+\dfrac{103}{96}\mathrm{e}^4\sin4a_m+\ldots
\]
\textit{La convergence de ce développement en série est assurée pour de faibles excentricités}(cf.Guy Sérane).

\section{L'animation à partir des données de la \textsc{Nasa}}
\begin{center}
\def\nFrames{230}% 230 images (229 jours)
\begin{animateinline}[controls,timeline=kepler16Nasa.dat,loop,%
                     begin={\begin{pspicture}(-8,-8)(8,8)},
                     end={\end{pspicture}}]{10}% 10 images/s
\pstVerb{%
/Pi 3.14159265359 def
/2Pi 6.28318530718 def
% les angles doivent être en degrés
% pour calculer les sinus, cosinus
% en postscript
/rad2deg { 180 mul Pi div } bind def % radians -> degrés
/deg2rad { 180 div Pi mul } bind def %  degrés -> radians
/tan {dup sin exch cos div} def
         /arctan {
            dup 0 ge
            {1 atan}
            {neg 1 atan neg}
         ifelse
        } def
/AV { % anomalie vraie
   1 exc add 1 exc sub div sqrt AE rad2deg 2 div tan mul arctan 2 mul
  } def
/mA 0.6897 def % masse étoile A
/mB 0.20255 def% masse étoile B
/mt mA mB add def % masse totale
/mP 3e-4 def % masse planète
/a1 0.22431 def % a pour le point réduit
/e1 0.15944 def % excentricité
/p1 a1 1 e1 dup mul sub mul def  % paramètre
/pA p1 mB mul mt div def % paramètre star A
/pB p1 mA mul mt div def % paramètre star B
/e2 0.0069 def % excentricité planète
/a2 0.7048 def % demi-grand axe planète
/w1 263.464 def % arguument
/w2 318 def % arguument planète
/p2 a2 1 e2 dup mul sub mul def  % paramètre planète
% t T e Kepler -> AE
% date periode excentricité -> Kepler -> anomalie excentrique
/Kepler {
  /exc exch def
  /Ti exch def
  /ti exch def
  /AM ti 2Pi mul Ti div def % anomalie moyenne à la date t
  /AE AM def % valeur initiale de anomalie excentrique
    5 {% 5 boucles
     /AETemp % anomalie excentrique provisoire
      AM exc rad2deg AE sin mul add
      def
     /AE AETemp def
    } repeat
    }def
% étoiles A et B
% T=41.076 jours e=0.15944
% planète
% T=228.776 j    e=0.0069
% on calcule sur la periode la planète
% on construit les tableaux des positions
% à 1 jour d'intervalle
/tabStarA [
 0 1 229 {/i exch def
  i 41.076 e1 Kepler
  /radius pA 1 e1 AV cos mul add div neg def
   radius AV w1 add cos mul
   radius AV w1 add sin mul
  } for
  ] def
/tabStarB [
 0 1 229 {/i exch def
  i 41.076 e1 Kepler
  /radius pB 1 e1 AV cos mul add div def
   radius AV w1 add cos mul
   radius AV w1 add sin mul
  } for
  ] def
 /tabPlanet [
 0 1 229 {/i exch def
  i 228.776 e2 Kepler
  /radius p2 1 e2 AV cos mul add div def
   radius AV w2 add cos mul
   radius AV w2 add sin mul
  } for
  ] def
    }%
\psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-8,-8)(8,8)
\psRandomStar[linecolor={[rgb]{1,1,0.5}},
    randomPoints=1000,sizeStar=1pt](-8,-8)(8,8){\psframe[linestyle=none](-8,-8)(8,8)}
%\listplot[unit=10,linecolor=gray]{tabStarA aload pop}
%\listplot[unit=10,linecolor=gray]{tabStarB aload pop}
%\listplot[unit=10,linecolor=gray]{tabPlanet aload pop}
\newframe
\multiframe{\nFrames}{i=0+2}{%
\pstVerb{/i1 \i\space def
  /xA tabStarA i1 get 10 mul def
  /yA tabStarA i1 1 add get 10 mul def
  /xB tabStarB i1 get 10 mul def
  /yB tabStarB i1 1 add get 10 mul def
  /xP tabPlanet i1 get 10 mul def
  /yP tabPlanet i1 1 add get 10 mul def
       }%
\pscircle[linestyle=none,fillstyle=gradient,gradmidpoint=0,gradend=yellow,GradientCircle=true,gradbegin=yellow!50!red!50](!xA yA){0.6}
\pscircle[linestyle=none,fillstyle=gradient,gradmidpoint=0,gradend=red,GradientCircle=true,gradbegin=red!50](!xB yB){0.2}
\pscircle*[linecolor=white](!xP yP){0.1}
}
\end{animateinline}
\end{center}
\end{document}

