Commit e143369d authored by Zoheyr Doctor's avatar Zoheyr Doctor 🔭
Browse files

Merge branch 'update_binary' into 'master'

Update binary

See merge request !1
parents 98969cf1 86924668
......@@ -46,7 +46,6 @@ def make_parser():
def main():
import json
import sys
import sys
# parse arguments
parser = make_parser()
args = parser.parse_args()
......@@ -97,8 +96,6 @@ def main():
binary = Binary(
m1 = m1,
m2 = m2,
c1 = 'k',
c2 = 'k',
alpha = 45.,
beta = 20.,
gamma = 80.,
......@@ -109,24 +106,8 @@ def main():
# evolve the binary over one orbit
n_frames = int(length*fps)
# populate positions for event
event_pos_dict = {
var: np.empty(n_frames)
for var in ['x1','y1','z1','x2','y2','z2']
}
binary.evolve(n_frames)
for i in range(n_frames):
binary.orbit()
pos1,pos2 = binary.project()
event_pos_dict['z1'][i] = pos1[0]
event_pos_dict['z2'][i] = pos2[0]
event_pos_dict['x1'][i] = pos1[1]
event_pos_dict['y1'][i] = pos1[2]
event_pos_dict['x2'][i] = pos2[1]
event_pos_dict['y2'][i] = pos2[2]
# write out evolution to dictionary
pos_dict[event] = event_pos_dict
binary_dict[event] = binary
artist_dict = {} # dictionary for matplotlib artists
......@@ -134,16 +115,19 @@ def main():
def init():
for ei,event in enumerate(events):
# create artists for both stars
# create artists for both stars
z1,x1,y1 = binary_dict[event].pos1_projected[0,:]
z2,x2,y2 = binary_dict[event].pos2_projected[0,:]
artist_dict['%s_m1'%event] = ax[ei].add_artist(
Circle((pos_dict[event]['x1'][0],pos_dict[event]['y1'][0]),
Circle((x1,y1),
radius=binary_dict[event].m1/200,
edgecolor='white',
facecolor='black'
)
)
artist_dict['%s_m2'%event] = ax[ei].add_artist(
Circle((pos_dict[event]['x2'][0],pos_dict[event]['y2'][0]),
Circle((x2,y2),
radius=binary_dict[event].m2/200,
edgecolor='white',
facecolor='black'
......@@ -151,7 +135,7 @@ def main():
)
# make sure star that is closer is shown on top of the other
zorder_1 = int(pos_dict[event]['z1'][0]>pos_dict[event]['z2'][0])
zorder_1 = int(z1>z2)
zorder_2 = ~zorder_1
artist_dict['%s_m1'%event].set_zorder(zorder_1)
artist_dict['%s_m2'%event].set_zorder(zorder_2)
......@@ -170,11 +154,13 @@ def main():
for ei,event in enumerate(events):
# update positions of stars
artist_dict['%s_m1'%event].set_center((pos_dict[event]['x1'][i],pos_dict[event]['y1'][i]))
artist_dict['%s_m2'%event].set_center((pos_dict[event]['x2'][i],pos_dict[event]['y2'][i]))
z1,x1,y1 = binary_dict[event].pos1_projected[i,:]
z2,x2,y2 = binary_dict[event].pos2_projected[i,:]
artist_dict['%s_m1'%event].set_center((x1,y1))
artist_dict['%s_m2'%event].set_center((x2,y2))
# make sure star that is closer is shown on top of the other
if pos_dict[event]['z1'][i]>pos_dict[event]['z2'][i]:
if z1>z2:
zorder_1 = 2
zorder_2 = 1
else:
......
#!/usr/bin/env python
import numpy as np
from twirl.binary import InspiralingBinary
from matplotlib import pyplot as plt
from matplotlib import animation
from matplotlib.patches import Circle
import sys
plt.style.use('dark_background')
fig,ax = plt.subplots()
fps = 30.
omega0 = 20.
m1 = 30.,
m2 = 30.,
binary = InspiralingBinary(
m1 = 30.,
m2 = 30.,
alpha = 45.,
beta = 20.,
gamma = 80.,
omega0 = omega0
)
print('tc:',binary.tc)
print('omegamax',binary.omegamax)
binary.inspiral()
print('t:',binary.t)
pos1 = binary.pos1_projected
pos2 = binary.pos2_projected
print('pos shape:',pos1.shape)
x1s = pos1[:,1]
y1s = pos1[:,2]
x2s = pos2[:,1]
y2s = pos2[:,2]
artist_dict = {}
def init():
artist_dict['m1'] = ax.add_artist(
Circle((x1s[0],y1s[0]),
radius=100.,
edgecolor='white',
facecolor='black'
)
)
artist_dict['m2'] = ax.add_artist(
Circle((x2s[0],y2s[0]),
radius=100.,
edgecolor='white',
facecolor='black'
)
)
ax.set_xlim([-2000,2000])
ax.set_ylim([-2000,2000])
return artist_dict
def animate(i):
artist_dict['m1'].set_center((x1s[i],y1s[i]))
artist_dict['m2'].set_center((x2s[i],y2s[i]))
return artist_dict
anim = animation.FuncAnimation(fig,animate,init_func=init,frames = len(x1s))
anim.save('inspiral.mp4',fps=fps,extra_args=['-vcodec', 'libx264'])
......@@ -73,12 +73,10 @@ def main():
) for event in events
}
pos_dict = {} # dictionary for all compact object positions
binary_dict = {} # dictionary for all Binary objects
# loop over all events
for event in events:
print('evolving ',event)
event_pos_dict = {} # dictionary for positions of this event
m1 = data[event]['mass_1_source']['med']
m2 = data[event]['mass_2_source']['med']
......@@ -86,8 +84,6 @@ def main():
binary = Binary(
m1 = m1,
m2 = m2,
c1 = 'k',
c2 = 'k',
alpha = 45.,
beta = 20.,
gamma = 80.,
......@@ -97,25 +93,9 @@ def main():
# evolve the binary over one orbit
n_frames = int(length*fps)
# populate positions for event
event_pos_dict = {
var: np.empty(n_frames)
for var in ['x1','y1','z1','x2','y2','z2']
}
for i in range(n_frames):
binary.orbit()
pos1,pos2 = binary.project()
event_pos_dict['z1'][i] = pos1[0]
event_pos_dict['z2'][i] = pos2[0]
event_pos_dict['x1'][i] = pos1[1]
event_pos_dict['y1'][i] = pos1[2]
event_pos_dict['x2'][i] = pos2[1]
event_pos_dict['y2'][i] = pos2[2]
binary.evolve(n_frames)
# write out evolution to dictionary
pos_dict[event] = event_pos_dict
binary_dict[event] = binary
......@@ -125,17 +105,20 @@ def main():
artist_dict = {} # dictionary for matplotlib artists
def init():
z1,x1,y1 = binary_dict[event].pos1_projected[0,:]
z2,x2,y2 = binary_dict[event].pos2_projected[0,:]
# create artists for both stars
artist_dict['%s_m1'%event] = ax.add_artist(
Circle((pos_dict[event]['x1'][0],pos_dict[event]['y1'][0]),
Circle((x1,y1),
radius=binary_dict[event].m1/200,
edgecolor='white',
facecolor='black'
)
)
artist_dict['%s_m2'%event] = ax.add_artist(
Circle((pos_dict[event]['x2'][0],pos_dict[event]['y2'][0]),
Circle((x2,y2),
radius=binary_dict[event].m2/200,
edgecolor='white',
facecolor='black'
......@@ -143,7 +126,7 @@ def main():
)
# make sure star that is closer is shown on top of the other
zorder_1 = int(pos_dict[event]['z1'][0]>pos_dict[event]['z2'][0])
zorder_1 = int(z1>z2)
zorder_2 = ~zorder_1
artist_dict['%s_m1'%event].set_zorder(zorder_1)
artist_dict['%s_m2'%event].set_zorder(zorder_2)
......@@ -155,13 +138,15 @@ def main():
return artist_dict
def animate(i):
# update positions of stars
artist_dict['%s_m1'%event].set_center((pos_dict[event]['x1'][i],pos_dict[event]['y1'][i]))
artist_dict['%s_m2'%event].set_center((pos_dict[event]['x2'][i],pos_dict[event]['y2'][i]))
z1,x1,y1 = binary_dict[event].pos1_projected[i,:]
z2,x2,y2 = binary_dict[event].pos2_projected[i,:]
artist_dict['%s_m1'%event].set_center((x1,y1))
artist_dict['%s_m2'%event].set_center((x2,y2))
# make sure star that is closer is shown on top of the other
if pos_dict[event]['z1'][i]>pos_dict[event]['z2'][i]:
if z1>z2:
zorder_1 = 2
zorder_2 = 1
else:
......
......@@ -13,6 +13,7 @@ setup(
scripts = [
'bin/animate_binaries_from_file',
'bin/make_event_gifs',
'bin/inspiral_test',
],
packages = [
'twirl',
......
......@@ -3,15 +3,9 @@ import numpy as np
class Binary:
def __init__(self,
m1=10.,
m2=10.,
c1='black',
c2='black',
alpha=50.,
beta=30.,
gamma=0.,
frame_rate = 30.,
omega = 1.
m1=10.,m2=10.,
alpha=50.,beta=30.,gamma=0.,
frame_rate = 30.,omega = 1.
):
"""
A class for representing binaries.
......@@ -20,8 +14,6 @@ class Binary:
----------
m1,m2: floats
masses of 1st and 2nd star, repsectively
c1,c2: matplotlib color
colors of 1st and 2nd star, respectively
inclination: float
inclination angle of binary wrt viewer
in degrees, between [0,90] deg.
......@@ -37,12 +29,10 @@ class Binary:
"""
self.m1 = m1
self.m2 = m2
self.c1 = c1
self.c2 = c2
self.alpha = alpha
self.beta = beta
self.gamma = gamma
self.omega = omega
self._omega = omega
self.frame_rate = frame_rate
self.precession='off'
self._phase1 = self.gamma_rad
......@@ -50,6 +40,18 @@ class Binary:
self._C = np.identity(3)
self._B = np.identity(3)
self._pos1_projected = []
self._pos2_projected = []
self._r = 1.
@property
def r(self):
return self._r
@property
def omega(self):
return self._omega
@property
def alpha_rad(self):
......@@ -108,15 +110,87 @@ class Binary:
"""
position of primary in XYZ coordinates of binary plane
"""
return ((self.m2/(self.m1+self.m2))
return (self.r*(self.m2/(self.m1+self.m2))
* np.array([np.cos(self.phase1),np.sin(self.phase1),0])
)
@property
def pos2(self):
return ((self.m1/(self.m1+self.m2))
return (self.r*(self.m1/(self.m1+self.m2))
* np.array([np.cos(self.phase2),np.sin(self.phase2),0])
)
def project(self):
return np.matmul(self.Amatrix,self.pos1),np.matmul(self.Amatrix,self.pos2)
return np.matmul(self.Amatrix,self.pos1),np.matmul(self.Amatrix,self.pos2)
def evolve(self,n_steps=1):
for i in range(n_steps):
self.orbit()
pos1_proj,pos2_proj = self.project()
self._pos1_projected.append(pos1_proj)
self._pos2_projected.append(pos2_proj)
@property
def pos1_projected(self):
return np.array(self._pos1_projected)
@property
def pos2_projected(self):
return np.array(self._pos2_projected)
class InspiralingBinary(Binary):
def __init__(self,
m1=10.,m2=10.,
alpha=50.,beta=30.,gamma=0.,
frame_rate = 30.,omega0 = 1.
):
super().__init__(
m1,m2,
alpha,beta,gamma,
frame_rate,omega0
)
self.Mc = ((self.m1*self.m2)**(3./5))/((self.m1+self.m2)**(1./5))
self.G = 4.3e-3 * 3.1e13 # km Msol^-1 (km/s)^2
self.c = 3e5 # km/s
self.tc = (((2.*omega0)**(-8./3))
* (5./((8.*np.pi)**(8./3)))
* (((self.c**3.)/(self.G*self.Mc))**(5./3))
)
print('tc',self.tc)
self.rmax = 2.*(self.G/(self.c**2))*(self.m1+self.m2)
self.omegamax = np.sqrt(self.G*(self.m1+self.m2)/self.rmax**3)
self.frame_rate = self.omegamax
self.t = [0.]
pos1_proj,pos2_proj = self.project()
self._pos1_projected = [pos1_proj]
self._pos2_projected = [pos2_proj]
@property
def omega(self):
return (0.5)*((((8.*np.pi)**(8./3))/5.)
* (self.G*self.Mc/(self.c**3.))**(5./3)
* (self.tc - self.t[-1])
)**(-3./8)
@property
def r(self):
return (self.G*(self.m1+self.m2)/(self.omega**2))**(1./3)
def inspiral(self):
while self.r > self.rmax:
print('r,rmax:',self.r,self.rmax)
self.t.append(self.t[-1]+(1./self.frame_rate))
self.orbit()
if isinstance(self.r,complex):
print('r going complex, breaking out of loop')
break
pos1_proj,pos2_proj = self.project()
self._pos1_projected.append(pos1_proj)
self._pos2_projected.append(pos2_proj)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment