Skip to content

Commit

Permalink
add debug statements, adjust SSB to LISA time conversions, fix sky co…
Browse files Browse the repository at this point in the history 
…ord labels
  • Loading branch information
@acorreia61201
acorreia61201 committed Sep 16, 2024
1 parent 533d253 commit 0e6c7d3
Showing 1 changed file with 52 additions and 27 deletions.
79 changes: 52 additions & 27 deletions pycbc/detector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -660,7 +660,7 @@ def overhead_antenna_pattern(right_ascension, declination, polarization):
from lisatools.detector import ESAOrbits
except ImportError:
ESAOrbits = None

try:
from fastlisaresponse import pyResponseTDI
except ImportError:
Expand Down Expand Up @@ -789,10 +789,10 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
The cross polarization of the GW in the radiation frame.
lamb : float
The ecliptic latitude of the source in the SSB frame.
The ecliptic longitude of the source in the SSB frame.
beta : float
The ecliptic longitude of the source in the SSB frame.
The ecliptic latitude of the source in the SSB frame.
polarization : float (optional)
The polarization angle of the GW in radians. Default 0.
Expand Down Expand Up @@ -820,7 +820,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,

# set waveform start time if specified
if reference_time is not None:
self.ref_time = reference_time
self.ref_time = reference_time + self.offset

# specify and cache sample times
start = self.ref_time
Expand All @@ -830,6 +830,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
hp.start_time = start
hc.start_time = start
self.sample_times = hp.sample_times.numpy()
print('cached wf inputs and times')

# make sure signal still lies within orbit length
assert hp.duration + start <= self.orbits.t_base[-1], (
Expand All @@ -842,14 +843,17 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,

# configure the orbit to match signal
self.orbits.configure(t_arr=self.sample_times)
print('configured orbits')

# rotate GW from radiation frame to SSB using polarization angle
hp, hc = self.apply_polarization(hp, hc, polarization)
print('applied polarization')

# format wf to hp + i*hc
hp = hp.numpy()
hc = hc.numpy()
wf = hp + 1j*hc
print('converted to numpy')

# save length of wf
self.n = len(wf)
Expand All @@ -868,14 +872,17 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,

if self.response_init is None:
# initialize the class
print('fresh init')
self.response_init = pyResponseTDI(1/self.dt, self.n, orbits=self.orbits,
use_gpu=use_gpu)
else:
# update params in the initialized class
print('update init')
self.response_init.sampling_frequency = 1/self.dt
self.response_init.num_pts = self.n
self.response_init.orbits = self.orbits
self.response_init.use_gpu = use_gpu
print('response initialized')

# project the signal
self.response_init.get_projections(wf, lamb, beta, t0=self.t0)
Expand Down Expand Up @@ -906,10 +913,10 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None,
The cross polarization of the GW in the radiation frame.
lamb : float
The ecliptic latitude in the SSB frame.
The ecliptic longitude in the SSB frame.
beta : float
The ecliptic longitude in the SSB frame.
The ecliptic latitude in the SSB frame.
polarization : float
The polarization angle of the GW in radians.
Expand Down Expand Up @@ -959,10 +966,16 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None,
# get dt from waveform data
self.dt = hp.delta_t

# save params for LISA time conversion
self.beta = beta
self.lamb = lamb
self.pol = polarization

# get index corresponding to time length t0
if t0 is not None:
self.t0 = t0
if pad_data or remove_garbage:
global pad_idx
pad_idx = int(self.t0/self.dt)

# pad the data with zeros
Expand All @@ -981,6 +994,7 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None,
# generate the Doppler time series
self.get_links(hp, hc, lamb, beta, polarization=polarization,
reference_time=reference_time, use_gpu=use_gpu)
print('get links')

# set TDI configuration (let FLR handle if not 1 or 2)
if tdi == 1:
Expand Down Expand Up @@ -1009,13 +1023,19 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None,

# generate the TDI channels
tdi_obs = self.response_init.get_tdi_delays()
print('tdi complete')

# processing
tdi_dict = {}
slc = slice(pad_idx, -pad_idx)
self.sample_times -= self.offset
self.ref_time -= self.offset
print('start preprocessing')

for i in range(len(tdi_chan)):
tdi_dict[tdi_chan[i]] = tdi_obs[i]
# save as TimeSeries with LISA frame times
tdi_dict[tdi_chan[i]] = TimeSeries(tdi_obs[i], delta_t=self.dt,
epoch=self.ref_time_LISA)
print(f'saved {i} to TimeSeries')

# treat start and end gaps
if remove_garbage:
Expand All @@ -1025,37 +1045,42 @@ def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None,
tdi_dict[tdi_chan[i]][-pad_idx:] = 0
elif type(remove_garbage) == bool:
# cut the edge data
slc = slice(pad_idx, -pad_idx)
tdi_dict[tdi_chan[i]] = tdi_dict[tdi_chan[i]][slc]
if i == 0:
# cut the sample times (only do once)
self.sample_times = self.sample_times[slc]
# update sample times once
self.sample_times = tdi_dict[tdi_chan[i]].sample_times
else:
raise ValueError('remove_garbage arg must be a bool ' +
raise ValueError('remove_garbage arg must be a bool ' +
'or "zero"')

# unapply offset for GPS SSB times
self.sample_times -= self.offset
self.ref_time -= self.offset

# save params for LISA time conversion
self.beta = beta
self.lamb = lamb
self.pol = polarization
print(f'finished postprocessing {i}')

return tdi_dict

@property
def sample_times_LISA(self):
def ref_time_LISA(self):
"""
Sample times of signal converted to LISA frame.
Signal start time converted to LISA frame.
"""
params = [self.sample_times, self.beta, self.lamb, self.pol]
params = [self.ref_time, self.lamb, self.beta, self.pol]
assert all(i is not None for i in params), ("Need to run project_wave for conversion")

# calculate first time converted to LISA; shift sample times to match LISA start time
lisa_start, _, _, _ = ssb_to_lisa(self.sample_times[0], self.beta, self.lamb,
self.pol, t0=0)
diff = lisa_start - self.sample_times[0]
# convert ref time to LISA
ssb_start = self.ref_time.gpsSeconds + 1e-9*self.ref_time.gpsNanoSeconds
lisa_start, _, _, _ = ssb_to_lisa(t_ssb = ssb_start,
longitude_ssb = self.lamb,
latitude_ssb = self.beta,
polarization_ssb = self.pol,
t0=self.offset)

return lisa_start

@property
def sample_times_LISA(self):
"""
Signal sample times converted to LISA frame.
"""
diff = self.ref_time_LISA - self.ref_time
return self.sample_times + diff


Expand Down

0 comments on commit 0e6c7d3

Please sign in to comment.