Skip to content

Commit

Permalink
long signals
Browse files Browse the repository at this point in the history
  • Loading branch information
@MarinerQ
MarinerQ committed Jul 11, 2024
1 parent b7a7aa7 commit e19369d
Show file tree
Hide file tree
Showing 4 changed files with 161 additions and 91 deletions.
3 changes: 2 additions & 1 deletion sealgw/calculation/localization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@
import lal
import lalsimulation
import bilby
#from gstlal import chirptime

# from gstlal import chirptime

# export OMP_NUM_THREADS=8
LAL_DET_MAP = dict(
Expand Down
78 changes: 39 additions & 39 deletions sealgw/simulation/antenna.py
View file
Original file line number Diff line number Diff line change
@@ -1,27 +1,25 @@
import numpy as np
import bilby
import bilby



class GWAntennaOnCPU():
class GWAntennaOnCPU:
def __init__(self, detector_tensor, vertex):
self.detector_tensor = detector_tensor
self.vertex = vertex

def getgha(self, gpstime, ra):
# Greenwich hour angle of source (radians).
gha = np.zeros_like(gpstime) - ra
for i,gpst in enumerate(gpstime):
for i, gpst in enumerate(gpstime):
gha[i] += bilby.gw.utils.greenwich_mean_sidereal_time(gpst)
return gha

def response(self, ra, dec, psi, gpstime):
bs = ra.shape[0]

X = np.zeros((bs, 3))
Y = np.zeros((bs, 3))


gha = self.getgha(gpstime, ra)

cosgha = np.cos(gha)
Expand All @@ -30,42 +28,44 @@ def response(self, ra, dec, psi, gpstime):
sindec = np.sin(dec)
cospsi = np.cos(psi)
sinpsi = np.sin(psi)

X[:,0] = -cospsi * singha - sinpsi * cosgha * sindec
X[:,1] = -cospsi * cosgha + sinpsi * singha * sindec
X[:,2] = sinpsi * cosdec

Y[:,0] = sinpsi * singha - cospsi * cosgha * sindec
Y[:,1] = sinpsi * cosgha + cospsi * singha * sindec
Y[:,2] = cospsi * cosdec

X[:, 0] = -cospsi * singha - sinpsi * cosgha * sindec
X[:, 1] = -cospsi * cosgha + sinpsi * singha * sindec
X[:, 2] = sinpsi * cosdec

Y[:, 0] = sinpsi * singha - cospsi * cosgha * sindec
Y[:, 1] = sinpsi * cosgha + cospsi * singha * sindec
Y[:, 2] = cospsi * cosdec

D = self.detector_tensor
fp = np.einsum('ij,jk,ik->i', X, D, X) - np.einsum('ij,jk,ik->i', Y, D, Y)
fc = np.einsum('ij,jk,ik->i', X, D, Y) + np.einsum('ij,jk,ik->i', Y, D, X)

return (fp, fc)

def time_delay_from_geocenter(self, ra, dec, gpstime):
bs = ra.shape[0]
gha = self.getgha(gpstime, ra)

cosgha = np.cos(gha)
singha = np.sin(gha)
cosdec = np.cos(dec)
sindec = np.sin(dec)

wavevector = np.zeros((bs, 3))
wavevector[:,0],wavevector[:,1],wavevector[:,2] = \
-cosgha*cosdec, cosdec*singha, -sindec

wavevector[:, 0], wavevector[:, 1], wavevector[:, 2] = (
-cosgha * cosdec,
cosdec * singha,
-sindec,
)

loc = self.vertex
dt = np.einsum('ij,j->i', wavevector, loc) / 299792458

return dt

def resp_and_dt(self, ra, dec, gpstime, psi):
bs = ra.shape[0]
bs = ra.shape[0]
X = np.zeros((bs, 3))
Y = np.zeros((bs, 3))

Expand All @@ -77,27 +77,27 @@ def resp_and_dt(self, ra, dec, gpstime, psi):
sindec = np.sin(dec)
cospsi = np.cos(psi)
sinpsi = np.sin(psi)

X[:,0] = -cospsi * singha - sinpsi * cosgha * sindec
X[:,1] = -cospsi * cosgha + sinpsi * singha * sindec
X[:,2] = sinpsi * cosdec

Y[:,0] = sinpsi * singha - cospsi * cosgha * sindec
Y[:,1] = sinpsi * cosgha + cospsi * singha * sindec
Y[:,2] = cospsi * cosdec
X[:, 0] = -cospsi * singha - sinpsi * cosgha * sindec
X[:, 1] = -cospsi * cosgha + sinpsi * singha * sindec
X[:, 2] = sinpsi * cosdec

Y[:, 0] = sinpsi * singha - cospsi * cosgha * sindec
Y[:, 1] = sinpsi * cosgha + cospsi * singha * sindec
Y[:, 2] = cospsi * cosdec



wavevector = np.zeros((bs, 3))
wavevector[:,0],wavevector[:,1],wavevector[:,2] = \
-cosgha*cosdec, cosdec*singha, -sindec


wavevector[:, 0], wavevector[:, 1], wavevector[:, 2] = (
-cosgha * cosdec,
cosdec * singha,
-sindec,
)

loc = self.vertex
D = self.detector_tensor

fp = np.einsum('ij,jk,ik->i', X, D, X) - np.einsum('ij,jk,ik->i', Y, D, Y)
fc = np.einsum('ij,jk,ik->i', X, D, Y) + np.einsum('ij,jk,ik->i', Y, D, X)
dt = np.einsum('ij,j->i', wavevector, loc) / 299792458
return (fp,fc,dt)

return (fp, fc, dt)
124 changes: 90 additions & 34 deletions sealgw/simulation/generating_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
)
from lal import LIGOTimeGPS
import logging
import lal
import lal

LAL_MTSUN_SI = lal.MTSUN_SI
LAL_PI = lal.PI
Expand Down Expand Up @@ -311,15 +311,15 @@ def f_of_tau(tau, m1=None, m2=None, mc=None):
return f


def tau_of_f( f, m1=None, m2=None, mc=None, chi=0):
def tau_of_f(f, m1=None, m2=None, mc=None, chi=0):
'''
Use 0PN if mc is provided. Otherwise use 3.5PN (TaylorF2), based on XLALSimInspiralTaylorF2ReducedSpinChirpTime
'''

if isinstance(f, list):
f = np.array(f)
if mc is None: # use 3.5PN

if mc is None: # use 3.5PN
if isinstance(f, (float, int, np.float64)):
tau = bilby.gw.utils.calculate_time_to_merger(f, m1, m2, safety=1)
elif isinstance(f, (np.ndarray)):
Expand All @@ -328,32 +328,61 @@ def tau_of_f( f, m1=None, m2=None, mc=None, chi=0):
eta = m1 * m2 / (m * m)
eta2 = eta * eta
chi2 = chi * chi
sigma0 = (-12769 * (-81. + 4. * eta)) / (16. * (-113. + 76. * eta) * (-113. + 76. * eta))
gamma0 = (565 * (-146597. + 135856. * eta + 17136. * eta2)) / (2268. * (-113. + 76. * eta))
sigma0 = (-12769 * (-81.0 + 4.0 * eta)) / (
16.0 * (-113.0 + 76.0 * eta) * (-113.0 + 76.0 * eta)
)
gamma0 = (565 * (-146597.0 + 135856.0 * eta + 17136.0 * eta2)) / (
2268.0 * (-113.0 + 76.0 * eta)
)

v = (LAL_PI * m * LAL_MTSUN_SI * f)**(1/3)
v = (LAL_PI * m * LAL_MTSUN_SI * f) ** (1 / 3)
tk = np.zeros((8, len(f))) # chirp time coefficients up to 3.5 PN

# chirp time coefficients up to 3.5PN
tk[0] = (5. * m * LAL_MTSUN_SI) / (256. * np.power(v, 8) * eta)
tk[1] = 0.
tk[2] = 2.9484126984126986 + (11 * eta) / 3.
tk[3] = (-32 * LAL_PI) / 5. + (226. * chi) / 15.
tk[4] = 6.020630590199042 - 2 * sigma0 * chi2 + (5429 * eta) / 504. + (617 * eta2) / 72.
tk[5] = (3 * gamma0 * chi) / 5. - (7729 * LAL_PI) / 252. + (13 * LAL_PI * eta) / 3.
tk[6] = -428.291776175525 + (128 * Pi_p2) / 3. + (6848 * LAL_GAMMA) / 105. + (3147553127 * eta) / 3.048192e6 - \
(451 * Pi_p2 * eta) / 12. - (15211 * eta2) / 1728. + (25565 * eta2 * eta) / 1296. + (6848 * np.log(4 * v)) / 105.
tk[7] = (-15419335 * LAL_PI) / 127008. - (75703 * LAL_PI * eta) / 756. + (14809 * LAL_PI * eta2) / 378.
tk[0] = (5.0 * m * LAL_MTSUN_SI) / (256.0 * np.power(v, 8) * eta)
tk[1] = 0.0
tk[2] = 2.9484126984126986 + (11 * eta) / 3.0
tk[3] = (-32 * LAL_PI) / 5.0 + (226.0 * chi) / 15.0
tk[4] = (
6.020630590199042
- 2 * sigma0 * chi2
+ (5429 * eta) / 504.0
+ (617 * eta2) / 72.0
)
tk[5] = (
(3 * gamma0 * chi) / 5.0
- (7729 * LAL_PI) / 252.0
+ (13 * LAL_PI * eta) / 3.0
)
tk[6] = (
-428.291776175525
+ (128 * Pi_p2) / 3.0
+ (6848 * LAL_GAMMA) / 105.0
+ (3147553127 * eta) / 3.048192e6
- (451 * Pi_p2 * eta) / 12.0
- (15211 * eta2) / 1728.0
+ (25565 * eta2 * eta) / 1296.0
+ (6848 * np.log(4 * v)) / 105.0
)
tk[7] = (
(-15419335 * LAL_PI) / 127008.0
- (75703 * LAL_PI * eta) / 756.0
+ (14809 * LAL_PI * eta2) / 378.0
)

vk = v.reshape((len(f), 1)) ** np.arange(8) # v^k
tau = (1+np.sum(tk[2:,:] * vk.T[2:,:], axis=0) ) * tk[0]
tau = (1 + np.sum(tk[2:, :] * vk.T[2:, :], axis=0)) * tk[0]
elif isinstance(m1, np.ndarray):
m = m1 + m2
eta = m1 * m2 / (m * m)
eta2 = eta * eta
chi2 = chi * chi
sigma0 = (-12769 * (-81. + 4. * eta)) / (16. * (-113. + 76. * eta) * (-113. + 76. * eta))
gamma0 = (565 * (-146597. + 135856. * eta + 17136. * eta2)) / (2268. * (-113. + 76. * eta))
sigma0 = (-12769 * (-81.0 + 4.0 * eta)) / (
16.0 * (-113.0 + 76.0 * eta) * (-113.0 + 76.0 * eta)
)
gamma0 = (565 * (-146597.0 + 135856.0 * eta + 17136.0 * eta2)) / (
2268.0 * (-113.0 + 76.0 * eta)
)

# Expand dimensions for batch processing
eta = eta[:, np.newaxis]
Expand All @@ -365,26 +394,53 @@ def tau_of_f( f, m1=None, m2=None, mc=None, chi=0):
m_expanded = m[:, np.newaxis]
f_expanded = f[np.newaxis, :]

v = (LAL_PI * m_expanded * LAL_MTSUN_SI * f_expanded)**(1/3)
tk = np.zeros((len(m1), 8, len(f))) # chirp time coefficients up to 3.5 PN for each batch
v = (LAL_PI * m_expanded * LAL_MTSUN_SI * f_expanded) ** (1 / 3)
tk = np.zeros(
(len(m1), 8, len(f))
) # chirp time coefficients up to 3.5 PN for each batch

# chirp time coefficients up to 3.5PN
tk[:, 0, :] = (5. * m_expanded * LAL_MTSUN_SI) / (256. * np.power(v, 8) * eta)
tk[:, 1, :] = 0.
tk[:, 2, :] = 2.9484126984126986 + (11 * eta) / 3.
tk[:, 3, :] = (-32 * LAL_PI) / 5. + (226. * chi) / 15.
tk[:, 4, :] = 6.020630590199042 - 2 * sigma0 * chi2 + (5429 * eta) / 504. + (617 * eta2) / 72.
tk[:, 5, :] = (3 * gamma0 * chi) / 5. - (7729 * LAL_PI) / 252. + (13 * LAL_PI * eta) / 3.
tk[:, 6, :] = -428.291776175525 + (128 * Pi_p2) / 3. + (6848 * LAL_GAMMA) / 105. + (3147553127 * eta) / 3.048192e6 - \
(451 * Pi_p2 * eta) / 12. - (15211 * eta2) / 1728. + (25565 * eta2 * eta) / 1296. + (6848 * np.log(4 * v)) / 105.
tk[:, 7, :] = (-15419335 * LAL_PI) / 127008. - (75703 * LAL_PI * eta) / 756. + (14809 * LAL_PI * eta2) / 378.

vk = np.power(v[:,:, np.newaxis], np.arange(8)).swapaxes(1,2) # v^k with correct broadcasting
tk[:, 0, :] = (5.0 * m_expanded * LAL_MTSUN_SI) / (
256.0 * np.power(v, 8) * eta
)
tk[:, 1, :] = 0.0
tk[:, 2, :] = 2.9484126984126986 + (11 * eta) / 3.0
tk[:, 3, :] = (-32 * LAL_PI) / 5.0 + (226.0 * chi) / 15.0
tk[:, 4, :] = (
6.020630590199042
- 2 * sigma0 * chi2
+ (5429 * eta) / 504.0
+ (617 * eta2) / 72.0
)
tk[:, 5, :] = (
(3 * gamma0 * chi) / 5.0
- (7729 * LAL_PI) / 252.0
+ (13 * LAL_PI * eta) / 3.0
)
tk[:, 6, :] = (
-428.291776175525
+ (128 * Pi_p2) / 3.0
+ (6848 * LAL_GAMMA) / 105.0
+ (3147553127 * eta) / 3.048192e6
- (451 * Pi_p2 * eta) / 12.0
- (15211 * eta2) / 1728.0
+ (25565 * eta2 * eta) / 1296.0
+ (6848 * np.log(4 * v)) / 105.0
)
tk[:, 7, :] = (
(-15419335 * LAL_PI) / 127008.0
- (75703 * LAL_PI * eta) / 756.0
+ (14809 * LAL_PI * eta2) / 378.0
)

vk = np.power(v[:, :, np.newaxis], np.arange(8)).swapaxes(
1, 2
) # v^k with correct broadcasting
tau = (1 + np.sum(tk[:, 2:, :] * vk[:, 2:, :], axis=1)) * tk[:, 0, :]
else:
print(f, type(f))
raise Exception("Unsupported type of f")
else: # use 0PN
else: # use 0PN
tau = 2.18 * (1.21 / mc) ** (5 / 3) * (100 / f) ** (8 / 3)
return tau

Expand Down
Loading

0 comments on commit e19369d

Please sign in to comment.