2 \documentclass{article
}
3 \usepackage[a4paper,margin=
2cm
]{geometry
}
4 \usepackage[T1]{fontenc}
5 \usepackage[latin1]{inputenc}%
6 \usepackage[garamond
]{mathdesign
}
7 \usepackage{pst-eqdf,pst-node,pst-tools
}
8 \usepackage{array,amsmath
}
11 \newpsstyle{vecteurA
}{arrowinset=
0.05,arrowsize=
0.1,linecolor=
{[rgb
]{1 0.5 0}}}
12 \newpsstyle{vecteurB
}{arrowinset=
0.05,arrowsize=
0.1,linecolor=
{[rgb
]{0 0.5 1}}}
13 \newpsstyle{vecteurC
}{arrowinset=
0.1,arrowsize=
0.2,linecolor=
{[rgb
]{1 0 0}}}
15 %% adapté de \psRandom du package pstricks-add
16 %% pour rendre aléatoire la taille des étoiles
19 \def\psset@sizeStar
#1{\pssetlength\pssizeStar{#1}}
21 \define@key
[psset
]{pst-eqd
}{randomPoints
}[1000]{\def\psk@randomPoints
{#1}}
22 \psset[pst-eqd
]{randomPoints=
1000}
23 \define@boolkey
[psset
]{pst-eqd
}[Pst@
]{color}[true
]{}
24 \psset[pst-eqd
]{color=false
}
25 \def\psRandomStar{\pst@object
{psRandomStar
}}%
26 \def\psRandomStar@i
{\@ifnextchar(
{\psRandomStar@ii
}{\psRandomStar@iii(
0,
0)(
1,
1)
}}
27 \def\psRandomStar@ii(
#1)
{\@ifnextchar(
{\psRandomStar@iii(
#1)
}{\psRandomStar@iii(
0,
0)(
#1)
}}
28 \def\psRandomStar@iii(
#1)(
#2)
#3{%
30 \ifx\pst@tempA
\pst@empty
\psclip{\psframe(
#2)
}\else\psclip{#3}\fi
31 \pst@getcoor
{#1}\pst@tempA
32 \pst@getcoor
{#2}\pst@tempB
35 \pst@tempA
\space /yMin exch def
37 \pst@tempB
\space /yMax exch def
41 rrand srand
% initializes the random generator
42 /getRandReal
{ rand
2147483647 div
} def
44 /DS
\pst@number
\pssizeStar\space getRandReal mul def
45 \@nameuse
{psds@
\psk@dotstyle
}
46 \ifPst@
color getRandReal
1 1 sethsbcolor
\fi
47 getRandReal dx mul xMin add
48 getRandReal dy mul yMin add
50 \ifx\psk@fillstyle
\psfs@solid fill
\fi stroke
60 \begin{filecontents
}{kepler16.dat
}
264 \title{Gravitation : une planète à deux soleils comme Tatooine celle de Luke Skywalker dans
\textit{La Guerre des étoiles
}}
265 \date{10 juillet
2\,
012}
268 \section{La représentation avec les données de la NASA
}
269 Ceci est une tentative pour essayer de schématiser une planète orbitant autour d'une étoile binaire, comme Kepler-
16b. On parle dans ce cas-là d'une planète
\textit{circumbinaire
}.
270 Le site de la NASA dédié à cette planète
\footnote{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. Voici les caractéristiques utiles pour le schéma et l'animation.
272 \begin{tabular
}{|ll|
}
274 \textbf{ Paramètres
}&
\textbf{Valeurs
}\\
\hline
275 \textit{Étoile A
}&\\
\hline
276 Masse, $M_A(M_
{\astrosun})$&
0.6897\\
277 Rayon, $R_A(R_
{\astrosun})$&
0.6489\\
\hline
278 \textit{Étoile B
}&\\
\hline
279 Masse, $M_B(M_
{\astrosun})$&
0.20255\\
280 Rayon, $R_B(R_
{\astrosun})$&
0.22623\\
\hline
281 \textit{Planète b
}&\\
\hline
282 Masse, $M_b(M_
{\jupiter})$&
0.333\\
283 Rayon, $R_b(R_
{\jupiter})$&
0.7538\\
\hline
284 \textit{Orbite de l'étoile binaire
}&\\
\hline
285 Période(jours)&
41.076\\
286 Demi-grand axe(ua) $
\mathrm{a_1
}$&
0.22431\\
287 Excentricité $
\mathrm{e_1
}$&
0.15944\\
288 Argument du périastre(deg) $
\omega_1$&
263.464\\
\hline
289 \textit{Orbite de la planète circumbinaire
} &\\
\hline
290 Période(jours)&
228.776\\
291 Demi-grand axe(ua)$
\mathrm{a_2
}$&
0.7048\\
292 Excentricité $
\mathrm{e_2
}$&
0.0069\\
293 Argument du périastre(deg) $
\omega_2$&
318\\
\hline
298 \begin{pspicture
}(-
5,-
5)(
5,
5)
299 \psgrid[subgriddiv=
0,gridcolor=lightgray,griddots=
10,gridlabels=
0pt
]%
300 \pstVerb{% pour le système binaire d'étoiles
301 /a1
0.22431 def
% semi-major axis
302 /e1
0.15944 def
% eccentricity
303 /w1
263.464 def
% argument of periaspe
304 /MA
0.6897 def
% masse of star A
305 /MB
0.20255 def
% masse of star B
306 /Mt MA MB add def
% masse totale
307 /P1
41.076 def
% period (day)
309 /Mb
0.3333 def
% masse of planet b (M_Jupiter)
310 /a2
0.7048 def
% semi-major axis
311 /e2
0.0069 def
% eccentricity
312 /w2
318 def
% argument of periaspe
313 /P2
228.776 def
% period (day)
314 % on en déduit le paramètre de chacune des ellipses
315 /p1 a1
1 e1 dup mul add mul def
316 /p2 a2
1 e2 dup mul add mul def
318 \pnode(!/radius p1
1 e1 w1 w1 sub cos mul add div def
319 /radius1 radius MB Mt div mul neg def
320 radius1 w1 cos mul
5 mul
321 radius1 w1 sin mul
5 mul)
{MA0
}
322 \parametricplot[plotpoints=
360,unit=
5,linecolor=red
]{0}{360}{
323 /radius p1
1 e1 t w1 sub cos mul add div def
324 /radius1 radius MB Mt div mul neg def
327 \pscircle[fillstyle=solid,fillcolor=yellow!
50](MA0)
{0.3}
328 \parametricplot[plotpoints=
360,unit=
5,linecolor=green
]{0}{360}{
329 /radius p1
1 e1 t w1 sub cos mul add div def
330 /radius2 radius MA Mt div mul def
333 \pnode(!/radius p1
1 e1 w1 w1 sub cos mul add div def
334 /radius1 radius MA Mt div mul def
335 radius1 w1 cos mul
5 mul
336 radius1 w1 sin mul
5 mul)
{MB0
}
337 \pscircle[fillstyle=solid,fillcolor=red
](MB0)
{0.15}
338 \parametricplot[plotpoints=
360,unit=
5,linecolor=blue
]{0}{360}{
339 /radius p2
1 e2 t w2 sub cos mul add div def
342 \pnode(!/radius p2
1 e2 w2 w2 sub cos mul add div def
343 radius w2 cos mul
5 mul
344 radius w2 sin mul
5 mul)
{Mb0
}
345 \pscircle[fillstyle=solid,fillcolor=blue
](Mb0)
{0.07}
348 \section{La simulation avec PSTricks
}
350 \begin{pspicture
}(-
3,-
1)(
6,
7)
351 \psgrid[subgriddiv=
0,gridcolor=lightgray,griddots=
10,gridlabels=
0pt
]%
352 \pstVerb{/x01 -
3 def /y01
0 def
353 /x02
5 def /y02
0 def
354 /xr012 x02 x01 sub def
355 /yr012 y02 y01 sub def
357 /xr002 x02 x0 sub def
358 /yr002 y02 y0 sub def
359 /xr001 x01 x0 sub def
360 /yr001 y01 y0 sub def
365 /xG012 x01 M1 mul x02 M2 mul add Mt div def
366 /yG012 y01 M1 mul y02 M2 mul add Mt div def
372 \pnode(!xG012 yG012)
{C12
} % centre de masse de M1 et M2
373 \pscircle[fillstyle=solid,fillcolor=yellow!
50](M2)
{0.21}
374 \pscircle[fillstyle=solid,fillcolor=red!
50](M1)
{0.5}
375 \pscircle[fillstyle=solid,fillcolor=blue!
50](M0)
{0.07}
376 \psline[linestyle=dotted
](M1)(M2)
377 \psline[linestyle=dashed
](M2)(M0)(M1)
378 \psline{<->
}(
6,
0)(
0,
0)(
0,
7)
382 \uput[l
](
0,
6.8)
{$
\mathcal{R
}$
}
383 \uput{0.6}[d
](M1)
{$M_1$
}
384 \uput{0.25}[d
](M2)
{$M_2$
}
385 \rput(M1)
{\psline[style=vecteurC
]{->
}(!xr012
3 div yr012
3 div)
\uput[d
](!xr012
3 div yr012
3 div)
{$
\overrightarrow{F
}_
{2/
1}$
}
386 \psline[style=vecteurC
]{->
}(!xr001
5 div neg yr001
5 div neg)
\uput[r
](!xr001
5 div neg yr001
5 div neg)
{$
\overrightarrow{F
}_
{0/
1}$
}}
387 \rput(M2)
{\psline[style=vecteurC
]{->
}(!xr012
3 div neg yr012
3 div neg)
\uput[d
](!xr012
3 div neg yr012
3 div neg)
{$
\overrightarrow{F
}_
{1/
2}$
}
388 \psline[style=vecteurC
]{->
}(!xr002
3 div neg yr002
3 div neg)
\uput[r
](!xr002
3 div neg yr002
3 div neg)
{$
\overrightarrow{F
}_
{0/
2}$
}}
389 \rput(M0)
{\psline[style=vecteurC
]{->
}(!xr002
3 div yr002
3 div)
\uput[r
](!xr002
3 div yr002
3 div)
{$
\overrightarrow{F
}_
{2/
0}$
}}
390 \rput(M0)
{\psline[style=vecteurC
]{->
}(!xr001
5 div yr001
5 div)
\uput[r
](!xr001
5 div yr001
5 div)
{$
\overrightarrow{F
}_
{1/
0}$
}}
393 Pour cela, on considère un système de trois corps en interaction gravitationnelle : l'étoile $M_1$ de masse $m_1$, l'étoile $M_2$ de masse $m_2$ constituant l'étoile binaire et une planète $P$ notée $M_0$, de masse $m$ orbitant autour de cette étoile double.
395 Avec : $
\overrightarrow{r_
{01}}=
\overrightarrow{M_0M_1
}$, $
\overrightarrow{r_
{02}}=
\overrightarrow{M_0M_2
}$ et $
\overrightarrow{r_
{12}}=
\overrightarrow{M_1M_2
}$ on pose :
397 \item $
\overrightarrow{r
}_
{01}=
\overrightarrow{r
}_1-
\overrightarrow{r
}_0$
398 \item $
\overrightarrow{r
}_
{02}=
\overrightarrow{r
}_2-
\overrightarrow{r
}_0$
399 \item $
\overrightarrow{r
}_
{12}=
\overrightarrow{r
}_2-
\overrightarrow{r
}_1$
401 Chacune des forces s'exprime par :
405 \overrightarrow{F
}_
{1/
0}=
\mathcal{G
}\dfrac{m_1m_0
}{r_
{10}^
3}\overrightarrow{r_
{01}}=-
\overrightarrow{F
}_
{0/
1}\\
[1em
]
406 \overrightarrow{F
}_
{2/
0}=
\mathcal{G
}\dfrac{m_2m_0
}{r_
{20}^
3}\overrightarrow{r_
{02}}=-
\overrightarrow{F
}_
{0/
2}\\
[1em
]
407 \overrightarrow{F
}_
{2/
1}=
\mathcal{G
}\dfrac{m_1m_2
}{r_
{12}^
3}\overrightarrow{r_
{12}}=-
\overrightarrow{F
}_
{1/
2}
411 L'application de la loi de Newton à chacun des corps donne :
415 m_0
\dfrac{\mathrm{d^
2}\overrightarrow{r_0
}}{\mathrm{d
}t^
2}=
\hphantom{-
}\mathcal{G
}\dfrac{m_1m_0
}{r_
{10}^
3}\overrightarrow{r_
{01}}+
416 \mathcal{G
}\dfrac{m_2m_0
}{r_
{20}^
3}\overrightarrow{r_
{02}}\\
[1em
]
417 m_1
\dfrac{\mathrm{d^
2}\overrightarrow{r_1
}}{\mathrm{d
}t^
2}=-
\mathcal{G
}\dfrac{m_0m_1
}{r_
{10}^
3}\overrightarrow{r_
{01}}+
418 \mathcal{G
}\dfrac{m_1m_2
}{r_
{12}^
3}\overrightarrow{r_
{12}}\\
[1em
]
419 m_2
\dfrac{\mathrm{d^
2}\overrightarrow{r_2
}}{\mathrm{d
}t^
2}=-
\mathcal{G
}\dfrac{m_0m_2
}{r_
{20}^
3}\overrightarrow{r_
{02}}-
420 \mathcal{G
}\dfrac{m_1m_2
}{r_
{12}^
3}\overrightarrow{r_
{12}}
424 Ce qui conduit à un système de
6 équations différentielles :
428 \ddot{x_0
}=
\hphantom{-
}\mathcal{G
}\dfrac{m_1
}{r_
{10}^
3}(x_1-x_0)+
\mathcal{G
}\dfrac{m_2
}{r_
{20}^
3}(x_2-x_0)\\
[1em
]
429 \ddot{y_0
}=
\hphantom{-
}\mathcal{G
}\dfrac{m_1
}{r_
{10}^
3}(y_1-y_0)+
\mathcal{G
}\dfrac{m_2
}{r_
{20}^
3}(y_2-y_0)\\
[1em
]
430 \ddot{x_1
}=-
\mathcal{G
}\dfrac{m_0
}{r_
{10}^
3}(x_1-x_0)+
\mathcal{G
}\dfrac{m_2
}{r_
{12}^
3}(x_2-x_1)\\
[1em
]
431 \ddot{y_1
}=-
\mathcal{G
}\dfrac{m_0
}{r_
{10}^
3}(y_1-y_0)+
\mathcal{G
}\dfrac{m_2
}{r_
{12}^
3}(y_2-y_1)\\
[1em
]
432 \ddot{x_2
}=-
\mathcal{G
}\dfrac{m_0
}{r_
{10}^
3}(x_2-x_0)-
\mathcal{G
}\dfrac{m_1
}{r_
{12}^
3}(x_2-x_1)\\
[1em
]
433 \ddot{y_2
}=-
\mathcal{G
}\dfrac{m_0
}{r_
{02}^
3}(y_2-y_0)-
\mathcal{G
}\dfrac{m_1
}{r_
{12}^
3}(y_2-y_1)
439 % 0 1 2 3 4 5 6 7 8 9 10 11
440 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
441 % x0 y0 x'0 y'0 x1 y1 x'1 y'1 x2 y2 x'2 y'2
442 \def\GravAlgIIIcorps{%
444 M1*(y
[4]-y
[0])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[8]-y
[0])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5|
%
445 M1*(y
[5]-y
[1])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[9]-y
[1])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5|
%
447 -M0*(y
[4]-y
[0])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[8]-y
[4])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
448 -M0*(y
[5]-y
[1])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[9]-y
[5])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
450 -M0*(y
[8]-y
[0])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5-M1*(y
[8]-y
[4])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
451 -M0*(y
[9]-y
[1])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5-M1*(y
[9]-y
[5])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5}
456 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
457 % x0 y0 x'0 y'0 x1 y1 x'1 y'1 x2 y2 x'2 y'2
458 \def\GravAlgIIIcorps{%
460 M1*(y
[4]-y
[0])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[8]-y
[0])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5|
%
461 M1*(y
[5]-y
[1])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[9]-y
[1])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5|
%
463 -M0*(y
[4]-y
[0])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[8]-y
[4])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
464 -M0*(y
[5]-y
[1])/((y
[4]-y
[0])^
2+(y
[5]-y
[1])^
2)^
1.5+M2*(y
[9]-y
[5])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
466 -M0*(y
[8]-y
[0])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5-M1*(y
[8]-y
[4])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5|
%
467 -M0*(y
[9]-y
[1])/((y
[8]-y
[0])^
2+(y
[9]-y
[1])^
2)^
1.5-M1*(y
[9]-y
[5])/((y
[8]-y
[4])^
2+(y
[9]-y
[5])^
2)^
1.5}
471 \def\InitCond{x00 y00 v0x0 v0y0 x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2
}
473 \begin{pspicture
}(-
8,-
8)(
8,
8)
474 \psframe*
[linecolor=
{[cmyk
]{1 1 0 0.7}}](-
8,-
8)(
8,
8)
478 dup mul neg
1 add sqrt
484 /x00
7 w2 cos mul def /y00
7 w2 sin mul def
493 /v0x0
0.36 w2 sin mul neg def
494 /v0y0
0.36 w2 cos mul def
499 /xG012 x01 M1 mul x02 M2 mul add Mt div def
500 /yG012 y01 M1 mul y02 M2 mul add Mt div def
502 \psgrid[subgriddiv=
0,gridcolor=white,griddots=
10,gridlabels=
0pt
](-
8,-
8)(
8,
8)
%
503 \psequadiff[plotpoints=
1000,algebraic,
515 saveData,filename=XAYA.dat
516 % ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
517 ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
518 \listplot[unit=
1,linecolor=yellow
]{XAYA aload pop
}
519 \psequadiff[plotpoints=
1000,algebraic,
531 saveData,filename=XBYB.dat
532 % ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
533 ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
534 \listplot[unit=
1,linecolor=red
]{XBYB aload pop
}
535 \psequadiff[plotpoints=
1000,algebraic,
547 saveData,filename=XPYP.dat
548 ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
549 \listplot[unit=
1,linecolor=blue
]{XbYb aload pop
}
550 \pnode(!x01 xG012 sub y01 yG012 sub)
{M01
}
551 \pnode(!x02 xG012 sub y02 yG012 sub)
{M02
}
552 \pnode(!x00 xG012 sub y00 yG012 sub)
{M00
}
553 \pscircle*
[linecolor=yellow
](M01)
{0.4}
554 \pscircle*
[linecolor=red
](M02)
{0.16}
555 \pscircle*
[linecolor=white
](M00)
{0.07}
556 %\rput(M01){\psline[unit=2,linecolor=red]{->}(!v0x1 v0y1)}
557 %\rput(M02){\psline[unit=2,linecolor=blue]{->}(!v0x2 v0y2)}
561 \section{Animation avec pst-eqdf et animate
}
563 \def\nFrames{200}% 200 images
564 \begin{animateinline
}[controls,timeline=kepler16.dat,loop,
%
565 begin=
{\begin{pspicture
}(-
8,-
8)(
8,
8)
},
566 end=
{\end{pspicture
}}]{5}% 5 images/s
567 \pstVerb{/XY1
[(XAYA.dat) run
] def
568 /XY2
[(XBYB.dat) run
] def
569 /XY3
[(XPYP.dat) run
] def
571 \psframe*
[linecolor=
{[cmyk
]{1 1 0 0.7}}](-
8,-
8)(
8,
8)
572 \psRandomStar[linecolor=
{[rgb
]{1,
1,
0.5}},
573 randomPoints=
1000,sizeStar=
1pt
](-
8,-
8)(
8,
8)
{\psframe[linestyle=none
](-
8,-
8)(
8,
8)
}
574 %\listplot[linecolor=gray]{XY1 aload pop}
575 %\listplot[linecolor=gray]{XY2 aload pop}
576 %\listplot[linecolor=gray]{XY3 aload pop}
578 \multiframe{\nFrames}{i=
0+
10}{% 1 point sur 10
579 \pstVerb{/X1 XY1
\i\space get def
580 /Y1 XY1
\i\space 1 add get def
581 /X2 XY2
\i\space get def
582 /Y2 XY2
\i\space 1 add get def
583 /X3 XY3
\i\space get def
584 /Y3 XY3
\i\space 1 add get def
587 \pscircle*
[linecolor=yellow
](!X1 Y1)
{0.5}
588 \pscircle*
[linecolor=red
](!X2 Y2)
{0.2}
589 \pscircle*
[linecolor=white
](!X3 Y3)
{0.07}