Skip to content

Commit

Permalink
edited polarization rotation function; final touches for first draft
Browse files Browse the repository at this point in the history
  • Loading branch information
@acorreia61201
acorreia61201 committed Jul 12, 2024
1 parent 2e21a1a commit 90daee5
Showing 1 changed file with 74 additions and 36 deletions.
110 changes: 74 additions & 36 deletions pycbc/detector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -648,7 +648,6 @@ def overhead_antenna_pattern(right_ascension, declination, polarization):


""" LISA class """

try:
import cupy
except ImportError:
Expand Down Expand Up @@ -694,7 +693,7 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use
# specify and cache the start time and orbital time series
if reference_time is None:
reference_time = self.orbits.t_base[0]
self.ref_time = reference_time
self.com_ref_time = reference_time
self.sample_times = None

# cache the FLR instance along with dt and n
Expand All @@ -707,25 +706,28 @@ def __init__(self, detector='LISA', reference_time=None, orbits=ESAOrbits(), use
# initialize whether to use gpu; FLR has handles if this cannot be done
self.use_gpu = use_gpu

def source_to_ssb(self, hp, hc, polarization):
def apply_polarization(self, hp, hc, polarization):
"""
Project a source-frame signal to SSB.
Apply polarization rotation matrix.
.. math::
\bmatrix{}
Parameters
----------
hp : pycbc.types.TimeSeries
The plus polarization of the GW in the radiation frame.
hp : array
The plus polarization of the GW.
hc : pycbc.types.TimeSeries
The cross polarization of the GW in the radiation frame.
hc : array
The cross polarization of the GW.
polarization : float
The polarization angle of the GW in radians.
Returns
-------
(pycbc.types.TimeSeries, pycbc.types.TimeSeries)
The plus and cross polarizations of the GW in the SSB frame.
(array, array)
The plus and cross polarizations of the GW rotated by the polarization angle.
"""
cphi = cos(2*polarization)
sphi = sin(2*polarization)
Expand All @@ -734,6 +736,44 @@ def source_to_ssb(self, hp, hc, polarization):
hc_ssb = hp*sphi + hc*cphi

return hp_ssb, hc_ssb

def time_delay_from_ssb(self, orbits, reference_time):
"""
Calculate the time delay from the SSB frame to the LISA COM frame.
Parameters
----------
orbits: lisatools.detector.Orbits
The orbital information of the satellites.
reference_time: float
The time in seconds in the SSB frame.
Returns
-------
float
The additive time delay factor between the SSB and LISA COM frame
at the given time.
"""
C_SI = constants.c.value

# configure orbits if not already
try:
orbits.t
except ValueError:
orbits.configure(linear_interp_setup=True)

# get positions of spacecraft at reference time
sc_x = []
for i in [1, 2, 3]:
sc_x.append(orbits.get_pos(reference_time, i))

# average sc positions to get COM position
com_vec = sum(sc_x)/len(sc_x)
com_dist = sum(com_vec*com_vec)**0.5

# time delay from SSB to LISA COM is distance/c
return com_dist/C_SI

def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
use_gpu=None, tdi=2):
Expand All @@ -757,9 +797,7 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
The ecleptic longitude of the source in the SSB frame.
polarization : float (optional)
The polarization angle of the GW in radians. If polarization is
0, the hp and hc inputs are treated as SSB frame projections.
Default 0.
The polarization angle of the GW in radians. Default 0.
reference_time : float (optional)
The GPS start time of the GW signal in the SSB frame. Default to
Expand All @@ -781,35 +819,35 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
"""
from fastlisaresponse import pyResponseTDI

# get dt and n from waveform data
dt = hp.delta_t # timestep (s)
n = len(hp) # number of points

# cache n, dt, pad length
# get dt from waveform data
dt = hp.delta_t
self.dt = dt
self.n = n
self.pad_idx = int(self.t0/self.dt)

# set waveform start time and cache time series
if reference_time is not None:
self.ref_time = reference_time
hp.start_time = self.ref_time
hc.start_time = self.ref_time
frame_delay = self.time_delay_from_ssb(self.orbits, reference_time)
self.com_ref_time = reference_time + frame_delay
hp.start_time = self.com_ref_time
hc.start_time = self.com_ref_time
self.sample_times = hp.sample_times.numpy()

# rescale the orbital time series to match waveform
### requires bugfix in LISAanalysistools to function
### awaiting response from dev team for formal release
self.orbits.configure(t_arr=self.sample_times)

# convert to SSB frame (if pol = 0, the signal is already in SSB frame)
if polarization != 0.:
hp, hc = self.source_to_ssb(hp, hc, polarization)
# rotate GW from radiation frame to SSB using polarization angle
hp, hc = self.apply_polarization(hp, hc, polarization)

# format wf to FLR standard hp + i*hc as ndarray
hp = hp.numpy()
hc = hc.numpy()
wf = hp + 1j*hc

# get n from waveform data
n = len(wf)
self.n = n

# set use_gpu to class input if not specified
if use_gpu is None:
Expand Down Expand Up @@ -853,8 +891,8 @@ def get_links(self, hp, hc, lamb, beta, polarization=0, reference_time=None,

return wf_proj

def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None, tdi=2,
tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=True):
def project_wave(self, hp, hc, lamb, beta, polarization, reference_time=None, tdi=2,
tdi_chan='AET', tdi_orbits=None, use_gpu=None, remove_garbage=False):
"""
Evaluate the TDI observables.
Expand All @@ -872,10 +910,8 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
beta : float
The ecleptic longitude in the SSB frame.
polarization : float (optional)
The polarization angle of the GW in radians. If polarization is
0, the hp and hc inputs are treated as SSB frame projections.
Default 0.
polarization : float
The polarization angle of the GW in radians.
reference_time : float (optional)
The GPS start time of the GW signal in the SSB frame. Default to
Expand Down Expand Up @@ -907,15 +943,16 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
dict
The TDI observables keyed by their corresponding TDI channel.
"""
start_time = time.time()

# set use_gpu
if use_gpu is None:
use_gpu = self.use_gpu

# generate the Doppler time series
self.get_links(hp, hc, lamb, beta, polarization=polarization, reference_time=reference_time,
use_gpu=use_gpu, tdi=tdi)

print(len(self.response_init.y_gw[0]))
print(len(self.sample_times))

# set TDI channels
if tdi_chan in ['XYZ', 'AET', 'AE']:
Expand Down Expand Up @@ -946,14 +983,15 @@ def project_wave(self, hp, hc, lamb, beta, polarization=0, reference_time=None,
if remove_garbage:
if remove_garbage == 'zero':
# zero the edge data
tdi_dict[tdi_chan[i]][self.pad_idx:] = 0
tdi_dict[tdi_chan[i]][:-self.pad_idx] = 0
tdi_dict[tdi_chan[i]][:self.pad_idx] = 0
tdi_dict[tdi_chan[i]][-self.pad_idx:] = 0
elif type(remove_garbage) == bool:
# cut the edge data
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]
print('check')
else:
# raise error
raise ValueError('remove_garbage arg must be a bool or "zero"')
Expand Down

0 comments on commit 90daee5

Please sign in to comment.