In [1]:
from __future__ import division, print_function
from ivisual import *
from math import *
In [13]:
scene=canvas(title="Binary Stars")
#constants
au=1.5e11
R=7e8
G=6.67e-11
#stars
m1=2e30
m2=2*m1
r1=vector(au,0,0)
r2=-m1/m2*r1
speed=1.1*3e4
v1=vector(0,speed,0)
v2=-m1/m2*v1
#time
t=0
dt=1e5
#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]
#ivisual objects
star1=sphere(radius=10*R, pos=r1, color=color.red, make_trail=True, retain=400)
star2=sphere(radius=10*R, pos=r2, color=color.yellow, make_trail=True, retain=400)
CM=sphere(radius=star1.radius/2, pos=(m1*r1+m2*r2)/(m1+m2), color=color.white)
while t<1000*dt:
rate(100)
r=r1-r2
rmag=mag(r)
runit=norm(r)
Fgrav=G*m1*m2/rmag/rmag
F1=Fgrav*(-runit)
v1=v1+F1/m1*dt
r1=r1+v1*dt
F2=-F1
v2=v2+F2/m2*dt
r2=r2+v2*dt
U=-G*m1*m2/rmag
K=1/2*m1*mag(v1)*mag(v1)+1/2*m2*mag(v2)*mag(v2)
E=U+K
Ulist.append(U)
Klist.append(K)
Elist.append(E)
rlist.append(rmag)
star1.pos=r1
star2.pos=r2
t=t+dt
In [5]:
import matplotlib.pyplot as plt
In [6]:
%matplotlib inline
In [7]:
#energy graphs
plt.title('energy vs. distance')
plt.xlabel('E (J)')
plt.ylabel('r (m)')
plt.plot(rlist,Ulist,'m.', rlist,Klist,'y.', rlist,Elist,'c.')
plt.show()
In this model, we will replace the two stars with a single star that has a mass $\mu$ that is equal to the reduced mass of the system, where
$$\mu=\frac{m_1m_2}{m_1+m_2}$$In this one-star model, its position is the relative position $\vec{r}$, where
$$\vec{r} = \vec{r}_1-\vec{r}_2$$$$r=|\vec{r}| = |\vec{r}_1-\vec{r}_2|$$Its radial velocity is $\dot{r}$.
Its angular momentum is
$$\vec{L}=\vec{r}_1 \times m_1\vec{v}_1+\vec{r}_2 \times m_2\vec{v}_2$$Its angular velocity $\dot{\phi}$ is
$$\dot{\phi}=\frac{|\vec{L}|}{r^2\mu}$$Its potential energy is
$$U=-\frac{Gm_1m_2}{r}+\frac{|\vec{L}|^2}{2\mu r^2}$$Its kinetic energy is
$$K=\frac{1}{2}\mu \dot{r}^2$$Its total energy is $E=K+U$. Note that because it is a single star in orbit (around nothing except a point in space!), it speeds up and slows down. Thus its kinetic energy increases and decreases and its potential energy increases and decreases. However its total energy $E$ remains constant.
The advantage of using a single star model is that it is completely described by two differential equations, one of which is first order.
$$\mu\ddot{r}=-\frac{Gm_1m_2}{r^2}+\frac{L^2}{\mu r^3}$$and
$$\mu \dot{\phi}=\frac{L}{r^2}$$where $L=|\vec{L}|$ is the magnitude of the angular momentum of the star.
The velocity of the single star is calculated in polar coordinates.
$$\vec{v}=<\dot{r},r\dot{\phi}>$$
In [10]:
scene2=canvas(title="One Star Model")
#reset variables
r1=vector(au,0,0)
r2=-m1/m2*r1
v1=vector(0,speed,0)
v2=-m1/m2*v1
#one star model
mu=m1*m2/(m1+m2)
L=m1*cross(r1,v1)+m2*cross(r2,v2)
Lmag=mag(L)
r=r1-r2
rmag=mag(r)
rdot=0
phi=0
phidot=Lmag/rmag/rmag/mu
v=vector(rdot,rmag*phidot)
#time
t=0
#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]
#ivisual objects
starmu=sphere(radius=20*R, pos=r, color=color.orange, make_trail=True, retain=400)
while t<1000*dt:
rate(100)
#update rmag
rdotdot=(-G*m1*m2/(rmag*rmag)+Lmag*Lmag/(mu*rmag*rmag*rmag))/mu
rdot=rdot+rdotdot*dt
rmag=rmag+rdot*dt
#update phi
phidot=Lmag/rmag/rmag/mu
phi=phi+phidot*dt
#calculate velocity and position
v=vector(rdot,rmag*phidot)
r=vector(rmag*cos(phi),rmag*sin(phi),0)
U=-G*m1*m2/rmag + Lmag*Lmag/(2*mu*rmag*rmag)
K=1/2*mu*rdot*rdot
E=U+K
Ulist.append(U)
Klist.append(K)
Elist.append(E)
rlist.append(rmag)
#ivisual objects
starmu.pos=r
t=t+dt
In [11]:
#energy graphs
plt.title('energy vs. distance')
plt.xlabel('E (J)')
plt.ylabel('r (m)')
plt.plot(rlist,Ulist,'m.', rlist,Klist,'y.', rlist,Elist,'c.')
plt.show()
Now here is the big deal! This single star and its motion can be used to tell you what the two-star system is doing without using Newton's second law.
Here is the same code as above. However, this time we will calculate the position of each star based on $r$ and $\phi$ and will draw them in the scene.
In [12]:
scene3=canvas(title="One Star Model with Two Stars Calculated")
#reset variables
r1=vector(au,0,0)
r2=-m1/m2*r1
v1=vector(0,speed,0)
v2=-m1/m2*v1
#one star model
mu=m1*m2/(m1+m2)
L=m1*cross(r1,v1)+m2*cross(r2,v2)
Lmag=mag(L)
r=r1-r2
rmag=mag(r)
rdot=0
phi=0
phidot=Lmag/rmag/rmag/mu
v=vector(rdot,rmag*phidot)
#time
t=0
#data lists
Ulist=[]
Klist=[]
Elist=[]
rlist=[]
#ivisual objects
starmu=sphere(radius=20*R, pos=r, color=color.orange, make_trail=True, retain=400)
star1=sphere(radius=10*R, pos=r1, color=color.red)
star2=sphere(radius=10*R, pos=r2, color=color.yellow)
CM=sphere(radius=star1.radius/2, pos=(m1*r1+m2*r2)/(m1+m2), color=color.white)
while t<1000*dt:
rate(100)
#update rmag
rdotdot=(-G*m1*m2/(rmag*rmag)+Lmag*Lmag/(mu*rmag*rmag*rmag))/mu
rdot=rdot+rdotdot*dt
rmag=rmag+rdot*dt
#update phi
phidot=Lmag/rmag/rmag/mu
phi=phi+phidot*dt
#calculate velocity and position
v=vector(rdot,rmag*phidot)
r=vector(rmag*cos(phi),rmag*sin(phi),0)
U=-G*m1*m2/rmag + Lmag*Lmag/(2*mu*rmag*rmag)
K=1/2*mu*rdot*rdot
E=U+K
Ulist.append(U)
Klist.append(K)
Elist.append(E)
rlist.append(rmag)
#ivisual objects
starmu.pos=r
star1.pos=m2/(m1+m2)*r
star2.pos=-m1/(m1+m2)*r
t=t+dt
In [ ]: