# ## 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
"""


