<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 BBH class
$ecc==0$
Step1:  5 realization results
(93328, 13)
(17287278, 13)
(17380606, 14)
Step2:  Mergers
Event Rate
Step3:  $\frac{1}{4\pi/3(30Mpc)^313.5Gyr}=\frac{1}{4\pi/32700013.5 {Gpc}^3 yr}$
Step4:  integrate from $z=0.6$,
which means t=7.7587 to 13.78, d_c=0 to 2.206 Gpc
cosmic calculation
Step5:  Consider cosmology, bhb merger in farther space can will be detectable
Step6:  Present BBHs
Step7:  $\dot{f}=k_0 f^{11/3}$
$k_0=\frac{96}{5}(2\pi)^{8/3}\frac{G^{5/3}}{c^5}\frac{m_1m_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/330^3 *10^{-9}$ Gpc^3 colume
```
Step8:  $f_N(19.6,10\sim100,10^{4.0})=[1190, 11897]$
in my case
Step9:  Generate file for LISA analysis
Parameters for binary source
$f_0$
Step10:  $T_{eject}$ variation for each model
Step11:  Weird models
48
210
<ASSISTANT_TASK:>
Python Code:
# ## 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]
dir_data='../data/Res/0220'
# 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))]
BBH['MC']=(BBH.M1*BBH.M2)**0.6*(BBH.M1+BBH.M2)**(-0.2)
Explanation: BBH class
$ecc==0$: never merge, period is lifetime in Gyr, seperation won't change
$ecc==1$: merge quickly, period is time of merge, seperation won't change
rest: 
merge: seperation==ecc==0, period is time of merge
present: period is in unit of second, seperation!=0, ecc!=0
End of explanation
display(merger.shape,present.shape,BBH.shape)
MC=(30*35)**0.6*(30+35)**(-0.2)
MC
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')
## Generate sample data for develop analysis scripts
# BBH.sample(10000).to_csv('../data/RES/1024/BBHex.dat')
Explanation: 5 realization results
(93328, 13)
(17287278, 13)
(17380606, 14)
End of explanation
merger.sample(5)
sum(merger.Ecc==1)
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')
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')
# sns.set_context("poster")
merge_curve=merge_kde.get_lines()[0].get_data()
Explanation: Mergers
Event Rate:
[LSC]BBH;aLIGO.pdf: $9 \sim 240Gpc^{-3}yr^{-1}$
[Carl;Morscher;Fred]BBH;GC;aLIGO.pdf: $10\sim100$ per year
My path:
merge event distribution over time
find the possibility of merge at recent, total merger: 93328 $(30Mpc)^{-3}(13.5Gyr)^{-1}$
scale it to $Gpc^{-3}yr^{-1}$
End of explanation
scale=len(merger)/(4.*pi/3.*30**3)/13.5
from scipy.interpolate import interp1d
f2=interp1d(merge_curve[0],merge_curve[1]*scale,kind='cubic')
f2(13.5)
plot(merge_curve[0],merge_curve[1]*scale)
title("merger event rate, per Gpc^3 per year, in local frame")
from scipy.integrate import quad
Explanation: $\frac{1}{4\pi/3(30Mpc)^313.5Gyr}=\frac{1}{4\pi/32700013.5 {Gpc}^3 yr}$
End of explanation
quad(f2, 7.7587,13.5)
Explanation: integrate from $z=0.6$,
which means t=7.7587 to 13.78, d_c=0 to 2.206 Gpc
cosmic calculation
End of explanation
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')
sns.jointplot(merger.Dist,merger.T_merge,kind='kde')
merger['MC']=(merger.M1*merger.M2)**0.6*(merger.M1+merger.M2)**(-0.2)
sns.distplot(merger.MC[merger.MC<50],kde=True)
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)
import matplotlib.pyplot as plt
import numpy as np
from moviepy.video.io.bindings import mplfig_to_npimage
import moviepy.editor as mpy
# 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)
import time
# 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)
Explanation: Consider cosmology, bhb merger in farther space can will be detectable
End of explanation
present.sample(5)
# 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)
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')
sum(log10(present.Period)<8.5)
present[log10(present.Period)<4]
present[log10(present.Period)<4.4]
present['MC']=(present.M1*present.M2)**0.6*(present.M1+present.M2)**(-0.2)
present.MC.describe()
# present[present.MC==min(present.MC)]
sns.distplot(present.MC)#[present.MC<50])
inspiral=present[log10(2./present.Period)>-4.0]
inspiral
sum(merger.Ecc==1)
merger.Seperation.describe()
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')
inspiral.ix[inspiral.index,['RA','Dec','Dist','Model','Age','T_eject','M1','M2','MC','Seperation','Ecc','Period']]
Explanation: Present BBHs
End of explanation
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
f_N(19.6,10,10**-4.0)
Explanation: $\dot{f}=k_0 f^{11/3}$
$k_0=\frac{96}{5}(2\pi)^{8/3}\frac{G^{5/3}}{c^5}\frac{m_1m_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/330^3 *10^{-9}$ Gpc^3 colume
```
End of explanation
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))
f_eta(19.6,10,10**-4.098)
Explanation: $f_N(19.6,10\sim100,10^{4.0})=[1190, 11897]$
in my case: $N=7$
$N=7$ only from dynamically formed bbh
field has $10^4$ more stars, might be [$10^{-2}$ to 0.5] less efficiency for bbh
bh formation merger events not counted
GC underestimated (GC SN \& dwarf Galaxy)
more bbh channel (galactic nucleus...)
End of explanation
present.shape
present.keys()
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
lisadata.sample(1000).describe().ix[['min','max']]
lisadata.to_csv('../data/RES/0220/lisa.dat.gz',
                sep='\t',index=None,compression='gzip')
lisadata[lisadata.Period<10**8.5].to_csv('../data/RES/0220/lisa_8.5.dat.gz',
                sep='\t',index=None,compression='gzip')
Explanation: Generate file for LISA analysis
Parameters for binary source
$f_0$: binary frequency
$\theta_s,\phi_s,\theta_L,\phi_L$: s for source; L for angular momentum vector
$\mathrm{A}=\frac{M_1M_2}{rD}$: amplitude
$\varphi_0$: phase at $t=0$
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
End of explanation
## 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)
# 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)  
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)
Explanation: $T_{eject}$ variation for each model
End of explanation
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()
bhblib[bhblib.Model.str.match(r"48-\d")]
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')
from gcpg import build
reload(build)
build.analysis('/Users/domi/Desktop/RES/1101')
Explanation: Weird models
48
210
End of explanation
<END_TASK>