Commit 9850db97 authored by Zoheyr Doctor's avatar Zoheyr Doctor 🔭
Browse files

made preliminary plot

parent 756155d9
......@@ -4,42 +4,114 @@ import numpy as np
from twirl.binary import Binary
from matplotlib import pyplot as plt
from matplotlib import animation
from matplotlib.patches import Circle
binary = Binary(
m1 = 10.,
m2 = 10.,
c1 = 'k',
c2 = 'k',
alpha = 45.,
beta = 20.,
gamma = 80.
)
x1 = []
y1 = []
x2 = []
y2 = []
for i in range(100):
binary.orbit()
pos1,pos2 = binary.project()
x1.append(pos1[1])
y1.append(pos1[2])
x2.append(pos2[1])
y2.append(pos2[2])
fig,ax = plt.subplots()
pt1, = ax.plot([],[],color='r',markersize=20,linestyle='none',marker='o')
def init():
pt1.set_data([],[])
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
return pt1,
def animate(i):
pt1.set_data([x1[i],x2[i]],[y1[i],y2[i]])
anim = animation.FuncAnimation(fig,animate,init_func=init,frames=100)
anim.save('orbit.mp4',fps=30,extra_args=['-vcodec', 'libx264'])
def make_parser():
from argparse import ArgumentParser
parser = ArgumentParser()
parser.add_argument(
"file",
type=str,
help='json file with median values',
)
return parser
def main():
import json
import sys
# parse arguments
parser = make_parser()
args = parser.parse_args()
with open(args.file,'r') as json_file:
data = json.load(json_file)
# get number of events and hence number of axes needed
events = list(data.keys())
n_events = len(events)
n_side = int(np.ceil(np.sqrt(n_events)))
n_side2 = int(np.ceil(n_events/float(n_side)))
fig,ax = plt.subplots(n_side,n_side2,figsize=(10,10))
ax = ax.flatten()
for axi in ax: axi.set_axis_off()
pos_dict = {}
binary_dict = {}
# loop over all events
for event in events:
print('evolving ',event)
event_pos_dict = {}
# create a Binary instance for the event
binary = Binary(
m1 = data[event]['mass_1_source']['med'],
m2 = data[event]['mass_2_source']['med'],
c1 = 'k',
c2 = 'k',
alpha = 45.,
beta = 20.,
gamma = 80.
)
# evolve the binary over one orbit
n_frames = 30
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]
# write out evolution to dictionary
pos_dict[event] = event_pos_dict
binary_dict[event] = binary
print(event_pos_dict['z1'],event_pos_dict['z2'])
artist_dict = {}
def init():
for ei,event in enumerate(events):
artist_dict['%s_m1'%event] = ax[ei].add_artist(Circle((event_pos_dict['x1'][0],event_pos_dict['y1'][0]),radius=binary_dict[event].m1/200,edgecolor='white',facecolor='black'))
artist_dict['%s_m2'%event] = ax[ei].add_artist(Circle((event_pos_dict['x2'][0],event_pos_dict['y2'][0]),radius=binary_dict[event].m2/200,edgecolor='white',facecolor='black'))
zorder_1 = int(event_pos_dict['z1'][0]>event_pos_dict['z2'][0])
zorder_2 = ~zorder_1
artist_dict['%s_m1'%event].set_zorder(zorder_1)
artist_dict['%s_m2'%event].set_zorder(zorder_2)
ax[ei].set_xlim([-1,1])
ax[ei].set_ylim([-1,1])
ax[ei].set_title(event)
return artist_dict
def animate(i):
for ei,event in enumerate(events):
artist_dict['%s_m1'%event].set_center((event_pos_dict['x1'][i],event_pos_dict['y1'][i]))
artist_dict['%s_m2'%event].set_center((event_pos_dict['x2'][i],event_pos_dict['y2'][i]))
if event_pos_dict['z1'][i]>event_pos_dict['z2'][i]:
zorder_1 = 2
zorder_2 = 1
else:
zorder_1 = 1
zorder_2 = 2
artist_dict['%s_m1'%event].set_zorder(zorder_1)
artist_dict['%s_m2'%event].set_zorder(zorder_2)
return artist_dict
anim = animation.FuncAnimation(fig,animate,init_func=init,frames=n_frames)
anim.save('orbit.gif',fps=30,dpi=300,writer='imagemagick',extra_args=['-vcodec', 'libx264'])
if __name__ == '__main__':
main()
......@@ -11,7 +11,7 @@ class Binary:
beta=30.,
gamma=0.,
frame_rate = 30.,
omega = 2.
omega = 1.
):
"""
A class for representing binaries.
......@@ -108,11 +108,15 @@ class Binary:
"""
position of primary in XYZ coordinates of binary plane
"""
return np.array([np.cos(self.phase1),np.sin(self.phase1),0])
return ((self.m2/(self.m1+self.m2))
* np.array([np.cos(self.phase1),np.sin(self.phase1),0])
)
@property
def pos2(self):
return np.array([np.cos(self.phase2),np.sin(self.phase2),0])
return ((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)
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