In [1]:
%pylab inline
In [2]:
import pandas as pd
from os import path
import seaborn as sns
sns.set_context("talk")
In [3]:
# ## load example data for testing
# BBHex=pd.read_csv('../data/RES/1024/BBHex.dat',delim_whitespace=True,header=None,
# names=['Galaxy','RA','Dec','Dist','VMag','Model','Age','T_eject','M1','M2','Seperation','Ecc','Period'])
# BBHex.sample(2)
# ## label example data for analysis
# merge_list=((BBHex.Seperation==0)&(BBHex.Ecc==0))|(BBHex.Ecc==1)
# present_list=-merge_list
# merger=BBHex[merge_list]
# merger=merger.rename(columns={'Period':'T_merge'})
# present=BBHex[present_list]
In [4]:
dir_data='../data/Res/0220'
In [5]:
# Load complete data set
BBH=pd.read_csv(path.join(dir_data,'BBHnow.dat.gz'),delim_whitespace=True,header=None,
names=['PGC','RA','Dec','Dist','VMag','Model',
'Age','T_eject','M1','M2','Seperation','Ecc','Period'])
merger=BBH[((BBH.Seperation==0)&(BBH.Ecc==0))|(BBH.Ecc==1)]
merger=merger.rename(columns={'Period':'T_merge'})
present=BBH[-(((BBH.Seperation==0)&(BBH.Ecc==0))|(BBH.Ecc==1))]
In [6]:
BBH['MC']=(BBH.M1*BBH.M2)**0.6*(BBH.M1+BBH.M2)**(-0.2)
5 realization results
(93328, 13)
(17287278, 13)
(17380606, 14)
In [7]:
display(merger.shape,present.shape,BBH.shape)
In [8]:
MC=(30*35)**0.6*(30+35)**(-0.2)
In [9]:
MC
Out[9]:
In [12]:
sns.distplot(BBH.MC[BBH.MC<50])
title('BBH merged: 93,328; present BBH: 17,287,278; total: 17,380,606')
xlabel('Chirp Mass')
annotate('GW150914',xy=(MC,0.018),xytext=(MC,0.04),arrowprops=dict(facecolor='red', shrink=0.05))
# savefig('../Fig/MC.jpeg',dpi=200,bbox_inches='tight')
Out[12]:
In [13]:
## Generate sample data for develop analysis scripts
# BBH.sample(10000).to_csv('../data/RES/1024/BBHex.dat')
In [14]:
merger.sample(5)
Out[14]:
In [15]:
sum(merger.Ecc==1)
Out[15]:
In [ ]:
In [16]:
sns.distplot(BBH.T_eject,kde=True,norm_hist=True)
text(7,0.3,BBH.T_eject.describe())
ylim([0,0.6])
xlim([-0.5,14])
# savefig(path.join(dir_data,'Fig/Teject_hist.jpeg'),dpi=200,bbox_inches='tight')
Out[16]:
In [17]:
merge_kde=sns.distplot(merger.T_merge,kde=True,norm_hist=True)
text(7,0.116,merger.T_merge.describe())
xlim([-0.5,14])
# savefig(path.join(dir_data,'Fig/Tmerge_hist.jpeg'),dpi=200,bbox_inches='tight')
Out[17]:
In [18]:
# sns.set_context("poster")
In [19]:
merge_curve=merge_kde.get_lines()[0].get_data()
$\frac{1}{4\pi/3*(30Mpc)^3*13.5Gyr}=\frac{1}{4\pi/3*27000*13.5 {Gpc}^3 yr}$
In [20]:
scale=len(merger)/(4.*pi/3.*30**3)/13.5
In [22]:
from scipy.interpolate import interp1d
f2=interp1d(merge_curve[0],merge_curve[1]*scale,kind='cubic')
f2(13.5)
Out[22]:
In [23]:
plot(merge_curve[0],merge_curve[1]*scale)
title("merger event rate, per Gpc^3 per year, in local frame")
Out[23]:
In [24]:
from scipy.integrate import quad
integrate from $z=0.6$, which means t=7.7587 to 13.78, d_c=0 to 2.206 Gpc
In [25]:
quad(f2, 7.7587,13.5)
Out[25]:
Consider cosmology, bhb merger in farther space can will be detectable
In [26]:
import numpy as np
import matplotlib.pyplot as pl
import scipy.stats as st
x = merger.Dist
y = merger.T_merge
xmin=5
xmax=35
ymin=0
ymax=15
# Peform the kernel density estimate
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([xx.ravel(), yy.ravel()])
values = np.vstack([x, y])
kernel = st.gaussian_kde(values)
f = np.reshape(kernel(positions).T, xx.shape)
fig = pl.figure()
ax = fig.gca()
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# Contourf plot
cfset = ax.contourf(xx, yy, f, cmap='Blues')
## Or kernel density estimate plot instead of the contourf plot
#ax.imshow(np.rot90(f), cmap='Blues', extent=[xmin, xmax, ymin, ymax])
# Contour plot
cset = ax.contour(xx, yy, f, colors='k')
# Label plot
ax.clabel(cset, inline=1, fontsize=10)
ax.set_xlabel('Distance [Mpc]')
ax.set_ylabel('Merge Time [Gyr]')
def tdelay(time,distance):
return time-distance*3.2625/1000
distance=np.arange(0,35,0.1)
y1=tdelay(13.5,distance)
y2=tdelay(13.4,distance)
# plt.fill_between(distance, y1, y2, color='grey', alpha='0.5')
# pl.savefig('../Fig/Event_rate.jpeg',dpi=200,bbox_inches='tight')
In [27]:
sns.jointplot(merger.Dist,merger.T_merge,kind='kde')
Out[27]:
In [28]:
merger['MC']=(merger.M1*merger.M2)**0.6*(merger.M1+merger.M2)**(-0.2)
In [29]:
sns.distplot(merger.MC[merger.MC<50],kde=True)
Out[29]:
In [31]:
fig = figure()#figsize=(12,8))
ax = fig.add_subplot(111, projection="mollweide")
# RA: [0h-24h] -> [-pi,pi]; Dec: [-90°,90°] -> [-pi/2,pi/2]\n"
ax.scatter((merger.RA-12)/12*pi,merger.Dec/180*pi,s=2)#,c=expGCpG.iMType,
# s=5*(2+log10(30./expGCpG.Dist.convert_objects(convert_numeric=True))),cmap=cm.Accent)
ax.set_xticklabels(['2h','4h','6h','8h','10h','12h','14h','16h','18h','20h','22h'])
ax.grid(True)
In [22]:
import matplotlib.pyplot as plt
import numpy as np
from moviepy.video.io.bindings import mplfig_to_npimage
import moviepy.editor as mpy
In [23]:
# duration = 1
# dt=0.2
# snap=lambda t:merger[(merger.T_merge>t)&(merger.T_merge<t+dt)]
# fig = figure()#figsize=(12,8))
# ax = fig.add_subplot(111, projection="mollweide")
# # RA: [0h-24h] -> [-pi,pi]; Dec: [-90°,90°] -> [-pi/2,pi/2]\n"
# ax.scatter((snap(0).RA-12)/12*pi,snap(0).Dec/180*pi)#,c=expGCpG.iMType,
# # s=5*(2+log10(30./expGCpG.Dist.convert_objects(convert_numeric=True))),cmap=cm.Accent)
# ax.set_xticklabels(['2h','4h','6h','8h','10h','12h','14h','16h','18h','20h','22h'])
# ax.grid(True)
# def make_frame_mpl(t):
# ax.scatter((snap(t/duration).RA-12)/12*pi,snap(t/duration).Dec/180*pi)
# return mplfig_to_npimage(fig) # RGB image of the figure
# animation =mpy.VideoClip(make_frame_mpl, duration=duration)
In [24]:
import time
In [25]:
# fig = figure()#figsize=(12,8))
# ax = fig.add_subplot(111, projection="mollweide")
# ax.set_xticklabels(['2h','4h','6h','8h','10h','12h','14h','16h','18h','20h','22h'])
# ax.grid(True)
# for t in linspace(0,13,260):
# # RA: [0h-24h] -> [-pi,pi]; Dec: [-90°,90°] -> [-pi/2,pi/2]\n"
# # ax.scatter((snap(0).RA-12)/12*pi,snap(0).Dec/180*pi)#,c=expGCpG.iMType,
# # s=5*(2+log10(30./expGCpG.Dist.convert_objects(convert_numeric=True))),cmap=cm.Accent)
# ax.scatter((snap(t).RA-12)/12*pi,snap(t).Dec/180*pi)
# time.sleep(0.01)
In [32]:
present.sample(5)
Out[32]:
In [33]:
# fig = figure()#figsize=(12,8))
# ax = fig.add_subplot(111, projection="mollweide")
# # RA: [0h-24h] -> [-pi,pi]; Dec: [-90°,90°] -> [-pi/2,pi/2]\n"
# ax.scatter((present.RA-12)/12*pi,present.Dec/180*pi,#c=expGCpG.iMType,
# s=log(2000000000/present.Period))#,cmap=cm.Accent)
# ax.set_xticklabels(['2h','4h','6h','8h','10h','12h','14h','16h','18h','20h','22h'])
# ax.grid(True)
In [34]:
sns.distplot(log10(2./(present.Period[log10(present.Period)<5.45])),rug=True,kde=True,bins=49)
xlim([-4.9,-3.5])
text(-4.2,1,(log10(2./(present.Period[log10(present.Period)<5.1])).describe()))
xlabel("Binary Frequency")
# savefig('../Fig/period.jpeg',dpi=200,bbox_inches='tight')
Out[34]:
In [121]:
sum(log10(present.Period)<8.5)
Out[121]:
In [36]:
present[log10(present.Period)<4]
Out[36]:
In [40]:
present[log10(present.Period)<4.4]
Out[40]:
In [41]:
present['MC']=(present.M1*present.M2)**0.6*(present.M1+present.M2)**(-0.2)
In [42]:
present.MC.describe()
Out[42]:
In [241]:
# present[present.MC==min(present.MC)]
In [355]:
sns.distplot(present.MC)#[present.MC<50])
Out[355]:
In [43]:
inspiral=present[log10(2./present.Period)>-4.0]
In [44]:
inspiral
Out[44]:
In [45]:
sum(merger.Ecc==1)
Out[45]:
In [46]:
merger.Seperation.describe()
Out[46]:
In [48]:
inspiral.ix[inspiral.index,['RA','Dec','Dist','Model','Age','T_eject','M1','M2','MC','Seperation','Ecc','Period']].to_csv('../data/RES/0220/present.dat',sep=' ',index=None,float_format='%.3f')
In [356]:
inspiral.ix[inspiral.index,['RA','Dec','Dist','Model','Age','T_eject','M1','M2','MC','Seperation','Ecc','Period']]
Out[356]:
$\dot{f}=k_0 f^{11/3}$
$k_0=\frac{96}{5}(2\pi)^{8/3}\frac{G^{5/3}}{c^5}\frac{m_1*m_2}{(m_1+m_2)^{1/3}}=\frac{96}{5}(2\pi)^{8/3}\frac{G^{5/3}}{c^5}M_{chirp}^{5/3}=const* M_{chirp}^{5/3}$
$const=3.68e-6 s^{5/3}$
$n=\frac{3\eta}{8k_0}f_{min}^{-8/3}$ where $f_{min}=\frac{2}{T_{max}\simeq 10^{4.4}[s]}\simeq 2*10^{-4.4} Hz$
N in [$\frac{4\pi}{3}30^3$ Mpc^3] = $\frac{3}{8k_0} f_{min}^{-8/3} \eta$ in [s x per Gpc^3 per year]
=> $N=\frac{3}{8k_0} f_{min}^{-8/3} \eta *[yr/s]*[Gpc^3/(V(30Mpc))]$
==> $\Delta t=\frac{3}{8k_0}f^{-8/3}_{now}$ per 30 Mpc^3 volume
when merge within 1 yr ($\Delta t \le 1yr$), $f \ge f_{min} = (\frac{8}{3} k_0 \Delta t)^{-3/8} $
$f_{min} = (\frac{8}{3}\frac{96}{5}(2\pi)^{8/3} \frac{G^{5/3}}{c^5}* 1yr)^{-3/8} M_{chirp}^{-5/8} \simeq 0.1164 M_{chirp}^{-5/8} Hz$ per 30 Mpc^3 volume, which is $4\pi/3*30^3 *10^{-9}$ Gpc^3 colume
In [49]:
def f_N(MC,eta,f_min):
return 3./8./(3.68e-6*MC**(5./3.))*(f_min)**(-8./3.)/3.154e7*(4.*3.14/3.*30.**3.*10**(-9))*eta
In [50]:
f_N(19.6,10,10**-4.0)
Out[50]:
in my case: $N=7$
In [51]:
def f_eta(MC,N,f_min):
return N*8./3.*(3.68e-6*MC**(5./3.))/(f_min)**(-8./3.)*3.154e7/(4.*3.14/3.*30.**3.*10**(-9))
In [52]:
f_eta(19.6,10,10**-4.098)
Out[52]:
In [117]:
Parameters for binary source
M1 = parameters(1)
M2 = parameters(2)
Porb = parameters(3) : orbital period, represent for orbital separation r, s
d = parameters(4) : distance between source and observer, pc
thetas = parameters(5) : based on galactic radec
phis = parameters(6) : based on galactic radec
thetal = parameters(7) : [0, pi] angular momentum
phil = parameters(8) : [0, 2pi] angular momentum
phio = parameters(9) : [0, 2pi] phase
gam = parameters(10) : [0, 2pi] amplitude modulation, Lambda
e = parameters(11) : eccentricity
In [53]:
present.shape
Out[53]:
In [54]:
present.keys()
Out[54]:
In [115]:
lisadata=present[['M1','M2','Period']] # M_sun, M_sun, second
lisadata['Dist']=present.Dist*1.0e6 # Mpc -> pc
lisadata['thetas']=present.Dec.apply(lambda x: (x+90.)/180.0*pi) # [-90,90]->[0,pi]
lisadata['phis']=present.RA.apply(lambda x: x*15.0/180.0*pi) # [0,24]->[0,2pi]
np.random.seed(seed=427)
lisadata['thetal']=np.random.rand(len(present))*pi # [0,pi]
lisadata['phil']=np.random.rand(len(present))*2*pi # [0,2pi]
lisadata['phio']=np.random.rand(len(present))*2*pi # [0,2pi]
lisadata['gam']=np.random.rand(len(present))*2*pi # [0,2pi]
lisadata['e']=present.Ecc
In [116]:
lisadata.sample(1000).describe().ix[['min','max']]
Out[116]:
In [118]:
lisadata.to_csv('../data/RES/0220/lisa.dat.gz',
sep='\t',index=None,compression='gzip')
In [124]:
lisadata[lisadata.Period<10**8.5].to_csv('../data/RES/0220/lisa_8.5.dat.gz',
sep='\t',index=None,compression='gzip')
In [ ]:
In [6]:
## Load BHLIB
bhlib=pd.DataFrame()
###################################
for i in range(1,325): # model
for j in range(1,11): # model id
###################################
bhe=pd.read_csv('../data/BHsystem/%d-%d-bhe.dat' %(i,j),
usecols=[0, 2, 3, 4, 6, 8, 10, 20], names=['T_eject','Type1','Type2','M1','M2','Seperation','Ecc','Model'],
header=None, delim_whitespace=True)
bhe.Model='%d-%d' %(i,j)
bhlib=pd.concat([bhlib,bhe],ignore_index=False)
In [7]:
# BHB binary
bhblib=bhlib[bhlib.Type1==14*(bhlib.Type2==14)].copy(deep=True)
bhsys=bhlib[-(bhlib.Type1==14*(bhlib.Type2==14))].copy(deep=True)
bhblib=bhblib.drop(bhblib.columns[[1, 2]], axis=1)
In [16]:
In [63]:
for i in range(1,325):
for j in range(1,6):
fig=figure(i)
try:
sns.distplot(bhblib.T_eject[bhblib.Model=="%d-%d" % (i,j)],rug=True,label='Model %d-%d' % (i,j))
except:
print "no ejected bbh in %d-%d" % (i,j)
legend()
savefig('../Fig/BHE_var/model-%d.jpeg' % i, dpi=200,bbox_inches='tight')
close(fig)
In [92]:
i=48
for j in range(1,6):
fig=figure(i)
try:
sns.distplot(bhblib.T_eject[bhblib.Model=="%d-%d" % (i,j)],kde=True,rug=True,label='Model %d-%d' % (i,j))
except:
print "no ejected bbh in %d-%d" % (i,j)
legend()
Out[92]:
In [104]:
bhblib[bhblib.Model.str.match(r"48-\d")]
Out[104]:
In [ ]:
In [ ]:
In [128]:
for i in range(1,6):
sns.distplot(bhblib.T_eject[bhblib.Model.str.contains('-%d' % i)],norm_hist=True,label='run %d' % i)
legend()
# savefig('../Fig/BHE_var/runs.jpeg',dpi=200,bbox_inches='tight')
Out[128]:
In [ ]:
In [23]:
from gcpg import build
In [26]:
reload(build)
Out[26]:
In [27]:
build.analysis('/Users/domi/Desktop/RES/1101')
In [ ]: