import numpy as np
from scipy.stats import truncnorm
from astropy import units as u
from .cosmology import Cosmology
from .constants import *
from .parameters import *
from .functions_earth import *
[docs]class BHB:
"""
This is a class to describe stellar mass black hole black hole mergers.
Parameters:
cosmos: define cosmological model
"""
def __init__(self, cosmos):
"""
Parameters:
cosmos (class): define cosmological model
"""
self.cosmos = cosmos
self.theta = thetaBHB
[docs] def set_model_theta(self, new_theta=None):
"""
The function to change theta from fiducial to user-defined in cosmic merger rate
Parameters:
new_theta (Optional[list of floats]): new parameter values
Returns:
(list of floats): theta with new values, if new_theta=None return theta to fiducial
"""
if new_theta != None:
self.theta = new_theta
else:
self.theta = thetaBHB
[docs] def mod_norm(self, Mch, F, iota, z):
"""
The normalization constant of the modulus of the GW waveform from stellar-mass binary black holes (BHB) mergers in frequency domain.
Parameters:
Mch (float): Red-shifted chirp mass of the BHB in solar masses
F (array of dtype float): The Antenna Patterns of the observatory
iota (float): The inclination of the BHB orbit
z (float): The redshift of the GW source
Returns:
(float): The normalization of the GW strain
"""
C = CONST_AMPLITUDE_GW*Mch**(5./6.)/self.cosmos.luminosity_distance(z).to_value(u.Gpc)
scaling_factor = np.sqrt(((1.+np.cos(iota)**2)/2.)**2*F[0]**2+np.cos(iota)**2*F[1]**2)
return C*scaling_factor
[docs] def mod_shape(self, f):
"""
The frequency dependence of the modulus of the GW waveform of BHB mergers.
Parameters:
f (float or array of dtype float): Frequencies of GW
Returns:
(float or array of dtype float): The modulus of the waveform corresponds to each frequency in f
"""
return f**(-7./6.)
[docs] def freq_limit(self, m1, m2, chi):
"""
The frequency upper limit of the waveform.
Parameters:
m1 (float): The Red-shifted individual masses
m2 (float): The Red-shifted individual masses
Returns:
(float): The upper limit of the frequency
"""
M = m1+m2
eta = sym_ratio(m1,m2)
y = np.zeros((4,3))
y[(1,0)] = 0.6437
y[(1,1)] = 0.827
y[(1,2)] = -0.2706
y[(2,0)] = -0.05822
y[(2,1)] = -3.935
y[(3,0)] = -7.092
sum_freq = 1.-4.455*(1-chi)**0.217+3.521**(1-chi)**0.26
for i in range(1,4):
N = min(3-i,2)
for j in range(0,N+1): ### N -> N+1
sum_freq += y[(i,j)]*eta**i*chi**j
return c**3/(GM_sun*M*np.pi)*sum_freq # unit [Hz]
[docs] def phase(self, f, m1, m2, chi, t0, phi0):
"""
The phase of the waveform of the stellar-mass BHB mergers in the frequency domain.
Parameters:
f (float or array of dtype float): Frequencies of GW
m1 (float): Red-shifted mass of the prime BH
m2 (float): Red-shifted mass of the second BH
chi (float): spin
t0 (float): The time of coalescence
phi0 (float): The phase of coalescence
Returns:
(array of dtype float): The phase psi(f) of the waveform as in exp^{−ipsi(f)}
"""
M = m1+m2
eta = sym_ratio(m1,m2)
x = np.zeros((4,3,8))
x[(1,0,2)] = -920.9
x[(1,1,2)] = 492.1
x[(1,2,2)] = 135.0
x[(2,0,2)] = 6742.
x[(2,1,2)] = -1053.
x[(3,0,2)] = -1.34e4
x[(1,0,3)] = 1.702e4
x[(1,1,3)] = -9566.
x[(1,2,3)] = 2182.
x[(2,0,3)] = -1.214e5
x[(2,1,3)] = 2.075e4
x[(3,0,3)] = 2.386e5
x[(1,0,4)] = -1.254e5
x[(1,1,4)] = 7.507e4
x[(1,2,4)] = 1.338e4
x[(2,0,4)] = 8.735e5
x[(2,1,4)] = -1.657e5
x[(3,0,4)] = -1.694e6
# x[(1,0,5)]=-920.9;x[(1,1,5)]=492.1;x[(1,2,5)]=135.0;x[(2,0,5)]=6742.;x[(2,1,5)]=-1053.;x[(3,0,5)]=-1.34e4
x[(1,0,6)] = -8.898e5
x[(1,1,6)] = 6.31e5
x[(1,2,6)] = 5.068e4
x[(2,0,6)] = 5.981e6
x[(2,1,6)] = -1.415e6
x[(3,0,6)] = -1.128e7
x[(1,0,7)] = 8.696e5
x[(1,1,7)] = -6.71e5
x[(1,2,7)] = -3.008e4
x[(2,0,7)] = -5.838e6
x[(2,1,7)] = 1.514e6
x[(3,0,7)] = 1.089e7
Psi = lambda chi : [0,0,3715./756.,-16.*np.pi+113.*chi*(1./3.),15293365.0/508032.0-405.*chi**2/8.0,0,0,0]
result = 1.
v = (np.pi*M*f*GM_sun/c**3)**(1./3.) # dimensional less, total mass
for k in range(2,8):
sum_psi = Psi(chi)[k]
for i in range(1,4):
N = min(3-i,2);
for j in range(0,N+1): ### N -> N+1
sum_psi += x[(i,j,k)]*eta**i*chi**j
result += sum_psi*v**k
return 2.*np.pi*f*t0+phi0+3./128./eta/v**5*result
[docs] def partial(self, A, m1, m2, chi, t0, phi0):
"""
The matrix of partial derivatives of waveform as function of f
Parameters:
A (float or array of dtype float): The modulus of h(f)
m1 (float): Red-shifted masses of the BHB
m2 (float): Red-shifted masses of the BHB
chi (float): spin
t0 (float): The time of coalescence
phi0 (float): The phase of coalescence
Returns:
(function): matrix of partial derivatives of (f,i,j), Pij(f)
"""
delta = 1e-5
phase_final = lambda f, m1, m2, chi, t0, phi0: self.phase(f, m1, m2, chi, t0, phi0)
# part_A=lambda f: 1/A;
# the error is determined with match-filtering of the waveform shape. therefore no derivative to the amplitude should be take.
part_m1 = lambda f:(self.phase(f, m1+delta*m1, m2, chi, t0, phi0)-self.phase(f, m1-delta*m1, m2, chi, t0, phi0))/(2.*delta*m1)
part_m2 = lambda f:(self.phase(f, m1, m2+delta*m2, chi, t0, phi0)-self.phase(f, m1, m2-delta*m2, chi, t0, phi0))/(2.*delta*m2)
part_chi = lambda f:(self.phase(f, m1, m2, chi+delta, t0, phi0)-self.phase(f, m1, m2, chi-delta, t0, phi0))/(2.*delta)
part_t0 = lambda f:2.*np.pi*f
part_phi = lambda f:1.
list_partial = [part_m1, part_m2, part_chi, part_t0, part_phi]
#PartMat = np.zeros((5,5))
#PartMat[(0,0)]=lambda f:1
#for i in range(1,5):
# for j in range(i,5):
# PartMat[(i,j)]=lambda f: list[i]*list[j]*A**2
return lambda f,i,j:list_partial[i](f)*list_partial[j](f)*self.mod_shape(f)**2*A**2
[docs] def cos_mer_rate(self, z, m1, m2, chi):
"""
The cosmic merger rate density of the stellar-mass BHB mergers.
Parameters:
z (float): The red-shift of the GW source
m1 (float): Red-shifted masses of the BHB
m2 (float): Red-shifted masses of the BHB
chi (float): spin
Returns:
(float): The number of mergers per year per Gpc3
"""
Mch = chirp_mass(m1,m2)
eta = sym_ratio(m1,m2)
a0,a1,t0,t1,b0,b1,c0,c1,d0,d1 = self.theta
if (z<=0 or z>=Z_HIGH or Mch<0 or eta<0 or eta>0.25):
result = 1e-10
norm = 1.
A = 1.
p_eta = 1.
p_chi = 1.
else:
t = self.cosmos.lookback_time(z).value #in unit of Gya
log10N = (a0+a1*t)/(np.exp((t-t0)/t1)+1.)
norm = 10**log10N
# the normalization
mu = b0+b1*t
sigL = c0+c1*t
sigR = d0+d1*t
A = np.sqrt(2./np.pi)/(sigL+sigR)
if Mch <= mu:
result = np.exp(-(mu-Mch)**2/(2.*sigL**2))
else:
result = np.exp(-(Mch-mu)**2/(2.*sigR**2))
# the distribution of chirpmass
q = (1/eta-2-np.sqrt((1/eta-2)**2-4))/2.
# the mass ratio m1/m2 (>1).
sig_q = 0.2
# the width of distribution of q
p_eta = np.exp(-(q-1)**2/(2.0*sig_q**2))
# the distribution of eta
sig_chi = 0.1
# the widht of distribution of chi
p_chi = np.exp(-chi**2/(2.0*sig_chi**2))
# the distribution of chi
return result*norm*A*p_eta*p_chi
# asymmetric (split) Gaussian function.
[docs] def tel_fun(self, z, m1, m2, chi, rho_cri, ant_fun, noise_tab):
"""
The telescope function of Laser Interferometers and kHz sources.
Parameters:
z (float): The redshift of the GW source
m1 (float): Red-shifted masses of the BHB
m2 (float): Red-shifted masses of the BHB
chi (float): spin
rho_cri (float): The detection SNR threshold
ant_fun (function): antenna pattern
noise_tab (array of dtype float): noise function for detector
Returns:
(float): The probability of detection
"""
# both masses should be intrinsic here.
length = 1000
f_up = self.freq_limit(m1=m1*(1+z), m2=m2*(1+z), chi=chi)
# the waveform want to eat red-shifted mass.
rho_sq_core_value = rho_sq_core(noise_tab, self.mod_shape, f_up=f_up)
theta_array = np.arccos(np.random.uniform(low=0., high=1., size=length))
varphi_array = np.random.uniform(low=0., high=2.*np.pi, size=length)
iota_array = np.arccos(np.random.uniform(low=0., high=1., size=length))
psi_array = np.random.uniform(low=0., high=2.*np.pi, size=length)
# an isotropic distribution of angles
F = ant_fun(theta_array, varphi_array, psi_array)
Mch = chirp_mass(m1,m2)
A_array = self.mod_norm(Mch*(1+z), F, iota_array, z)
rho_sq_array = np.array(4.*A_array**2*rho_sq_core_value)
heav_array = np.heaviside(rho_sq_array-rho_cri**2,0)
return np.mean(heav_array)
#return rho_sq_array
[docs] def density_det(self, T, z, m1, m2, chi, rho_cri, ant_fun, noise_tab):
"""
The parameter distribution of detected sources.
Parameters:
T (float): Observation time, in unit of minute
z (float): The redshift of the GW source
m1 (float): Red-shifted masses of the BHB
m2 (float): Red-shifted masses of the BHB
chi (float): spin
rho_cri (float): The detection SNR threshold
ant_fun (function): antenna pattern
noise_tab (array of dtype float): noise function for detector
Returns:
(float): Number density of detection
"""
T_year = (T*u.min).to(u.a).value # T in unit of minuts
return T_year*self.cos_mer_rate(z, m1, m2, chi)*self.tel_fun(z,m1,m2,chi, rho_cri,ant_fun,noise_tab)*self.cosmos.differential_comoving_volume(z).to_value(u.Gpc**3/u.sr)
[docs] def tot_num(self, T, rho_cri, ant_fun, noise_tab, ranges):
"""
The total number of sources detected.
Parameters:
T (float): Observation time, in unit of minutes
rho_cri (float): The detection SNR threshold
ant_fun (function): antenna pattern
noise_tab (array of dtype float): noise function for detector
ranges (list): parameter ranges
Returns:
(float): The total number of sources detected
"""
length = 5000
z_ranges = ranges[0]
m1_ranges = ranges[1]
m2_ranges = ranges[2]
chi_ranges = ranges[3]
z_array = np.random.uniform(low=z_ranges[0], high=z_ranges[1], size=length)
m1_array = np.random.uniform(low=m1_ranges[0], high=m1_ranges[1], size=length)
m2_array = np.random.uniform(low=m2_ranges[0], high=m2_ranges[1], size=length)
chi_array = np.random.uniform(low=chi_ranges[0], high=chi_ranges[1], size=length)
Ndet_array = np.array([self.density_det(T, z, m1, m2, chi,rho_cri,ant_fun,noise_tab) for (z, m1, m2, chi) in zip(z_array, m1_array, m2_array, chi_array)])
z_scope = z_ranges[1]-z_ranges[0]
m1_scope = m1_ranges[1]-m1_ranges[0]
m2_scope = m2_ranges[1]-m2_ranges[0]
chi_scope = chi_ranges[1]-chi_ranges[0]
return z_scope*m1_scope*m2_scope*chi_scope*np.mean(Ndet_array)
[docs] def MCMCsample(self, n_sources, inits, ranges, T, rho_cri, ant_fun, noise_tab):
"""
The Markov-chain-Monte-Carlo sampling according to the distribution ndet(M, eta, z).
Parameters:
n_sources (int): number of sources to detect
inits (list of floats): initial values of z, m1, m2, chi
ranges (list of list of floats): parameter ranges
rho_cri (float): The detection SNR threshold
ant_fun (function): antenna pattern
noise_tab (array of dtype float): noise function for detector
Returns:
(list of arrays of dtype float): samples of parameters Mchirp,z,m1,m2,chi,D
"""
burning_steps = 1000
jump_steps = 10
z0,m10,m20,chi0 = inits
step_z = 0.01
step_m1 = 1.
step_m2 = 1.
step_chi = 0.1
z_ranges = ranges[0]
m1_ranges = ranges[1]
m2_ranges = ranges[2]
chi_ranges = ranges[3]
# MCMC sampling according to distribution function(z,m1,m2)
z,m1,m2,chi = [z0,m10,m20,chi0]
# inital parameters
length = 1
samples_z = np.array([z])
samples_m1 = np.array([m1])
samples_m2 = np.array([m2])
samples_chi = np.array([chi])
while length < (n_sources*jump_steps+burning_steps):
z_next,m1_next,m2_next,chi_next = [z,m1,m2,chi]+np.random.normal(scale=[step_z,step_m1,step_m2,step_chi],size=4)
if z_ranges[0]<z_next<z_ranges[1] and m1_ranges[0]<m1_next<m1_ranges[1] and m2_ranges[0]<m2_next<m2_ranges[1] and chi_ranges[0]<chi_next<chi_ranges[-1]:
p0 = self.density_det(T,z,m1,m2,chi,rho_cri,ant_fun,noise_tab)
pac = min(1,self.density_det(T,z_next,m1_next,m2_next,chi_next,rho_cri,ant_fun,noise_tab)/p0) # accept probability
if pac==1 or np.random.uniform(low=0.0, high=1.0, size=None) <= pac:
z,m1,m2,chi = [z_next,m1_next,m2_next,chi_next]
samples_z = np.append(samples_z,[z])
samples_m1 = np.append(samples_m1,[max(m1,m2)])
samples_m2 = np.append(samples_m2,[min(m1,m2)])
samples_chi = np.append(samples_chi,[chi])
length += 1
samples_D = self.cosmos.luminosity_distance(samples_z).value
samples_Mc = chirp_mass(samples_m1,samples_m2)
return [samples_Mc[burning_steps::jump_steps],samples_z[burning_steps::jump_steps],samples_m1[burning_steps::jump_steps],samples_m2[burning_steps::jump_steps],samples_chi[burning_steps::jump_steps],samples_D[burning_steps::jump_steps]]
[docs] def errors_FIM(self, n_sources, samples, noise_tab):
"""
Return errors of parameters from Fisher Information Matrix
Parameters:
n_sources (int): number of sources to detect
samples (list of arrays of dtype float): sampled parameters
noise_tab (array of dtype float): noise function for detector
Returns:
(list of arrays of dtype float): errors of sampled parameters m1,m2,chi
"""
errors_m1 = np.zeros(n_sources)
errors_m2 = np.zeros(n_sources)
errors_chi = np.zeros(n_sources)
for i in range(0,n_sources):
z = samples[1][i]
m1 = samples[2][i]
m2 = samples[3][i]
Mch = samples[0][i]
DL = self.cosmos.luminosity_distance(z).to_value(u.Gpc)
A = CONST_AMPLITUDE_GW*(Mch*(1.+z))**(5./6.)/DL
chi = samples[4][i]
part_mat = self.partial(A,(1.+z)*m1,(1.+z)*m2,chi,0,0)
f_up = self.freq_limit((1.+z)*m1,(1.+z)*m2,chi)
fim = fis_inf_matr(part_mat,noise_tab,f_up=f_up,numpar=5)
err_mt = errors(fim)
#ErrMt_m1m2=ErrMt[np.ix_([0,1],[0,1])]
conv = FIM_conv_matr(m1,m2,5)
err_mt_tilde = np.dot(np.dot(conv.transpose(),err_mt),conv)
errors_m1[i] = np.sqrt(np.diag(err_mt))[0]
errors_m2[i] = np.sqrt(np.diag(err_mt))[1]
errors_chi[i] = np.sqrt(np.diag(err_mt))[2]
return [errors_m1,errors_m2,errors_chi]
[docs]class DNS(BHB):
"""
This is a class to describe double neutron stars mergers.
It inherits all functions from BHB class except set_model_theta and cos_mer_rate.
Attributes:
cosmos: define cosmological model
"""
def __init__(self,cosmos):
"""
Parameters:
cosmos (class): define cosmological model
"""
self.cosmos = cosmos
self.theta = thetaDNS
[docs] def set_model_theta(self, new_theta=None):
"""
The function to change theta from fiducial to user-defined.
Parameters:
new_theta (Optional[list of floats]): new parameter values
Returns:
(list of floats): theta with new values, if new_theta=None return theta to fiducial
"""
if new_theta != None:
self.theta = new_theta
else:
self.theta = thetaDNS
[docs] def cos_mer_rate(self, z, m1, m2, chi):
"""
The cosmic merger rate density of the stellar-mass DNS mergers.
Parameters:
z (float): The redshift of the GW source
m1 (float): Red-shifted masses of the NS
m2 (float): Red-shifted masses of the NS
chi (float): spin
Returns:
(float): The number of mergers per year per Gpc3
"""
a0,a1,t0,t1,loc,scal,low,high = self.theta
if (z<=0 or z>=Z_HIGH): rate = 1e-10
else:
t = self.cosmos.lookback_time(z).value
log10N = (a0+a1*t)/(np.exp((t-t0)/t1)+1.)
norm = 10**log10N
# the normalization
aNS,bNS = (low-loc)/scal, (high-loc)/scal
pm1 = truncnorm.pdf(m1, aNS, bNS, loc=loc, scale=scal)
pm2 = truncnorm.pdf(m2, aNS, bNS, loc=loc, scale=scal)
# the distribution of m1 and m2.
rate = norm*pm1*pm2
return rate
[docs]class BHNS(BHB):
"""
This is a class to describe stellar mass black hole neutron stars mergers.
It inherits all functions from BHB class except set_model_theta and cos_mer_rate.
Attributes:
cosmos: define cosmological model
"""
def __init__(self, cosmos):
"""
Parameters:
cosmos (class): define cosmological model
"""
self.cosmos = cosmos
self.theta = thetaBHNS
[docs] def set_model_theta(self, new_theta=None):
"""
The function to change theta from fiducial to user-defined.
Parameters:
new_theta (Optional[list of floats]): new parameter values
Returns:
(list of floats): theta with new values, if new_theta=None return theta to fiducial
"""
if new_theta != None:
self.theta = new_theta
else:
self.theta = thetaBHNS
[docs] def cos_mer_rate(self, z, M, m, chi):
"""
The cosmic merger rate density of the stellar-mass BHB mergers.
Parameters:
z (float): The redshift of the GW source
M (float): Red-shifted mass of the BH
m (float): Red-shifted mass of the NS
chi (float): spin
Returns:
(float): The number of mergers per year per Gpc3
"""
a0,a1,t0,t1,loc,scal,low,high,b0,b1,c0,c1,d0,d1 = self.theta
if (z<=0 or z>=Z_HIGH or M<0 or m<0): rate = 1e-10
else:
t = self.cosmos.lookback_time(z).value
log10N = (a0+a1*t)/(np.exp((t-t0)/t1)+1.)
norm = 10**log10N
# the normalization
aNS,bNS = (low-loc)/scal, (high-loc)/scal
pm = truncnorm.pdf(m, aNS, bNS, loc=loc, scale=scal)
# the mass of neutron star
mu = b0+b1*t
sigL = c0+c1*t
sigR = d0+d1*t
A = np.sqrt(2./np.pi)/(sigL+sigR)
if M <= mu: pM = np.exp(-(mu-M)**2/(2.*sigL**2))
else: pM = np.exp(-(M-mu)**2/(2.*sigR**2))
# the distribution of BH mass
rate = norm*pm*A*pM
return rate