Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use SkyGrid functions to write the HDF5 file and fix error in angular convention #5000

Merged
merged 9 commits into from
Jan 9, 2025
35 changes: 13 additions & 22 deletions bin/pycbc_make_sky_grid
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,8 @@ from scipy.spatial.transform import Rotation as R

import pycbc
from pycbc.detector import Detector
from pycbc.io.hdf import HFile
from pycbc.types import angle_as_radians
from pycbc.tmpltbank.sky_grid import SkyGrid


def spher_to_cart(sky_points):
Expand All @@ -36,6 +36,7 @@ def cart_to_spher(sky_points):
spher = np.zeros((len(sky_points), 2))
spher[:, 0] = np.arctan2(sky_points[:, 1], sky_points[:, 0])
spher[:, 1] = np.arcsin(sky_points[:, 2])
spher[spher[:, 0] <0, 0] += 2 * np.pi
return spher


Expand Down Expand Up @@ -163,25 +164,15 @@ rota = r.apply(cart)
# Convert cartesian coordinates back to spherical coordinates
spher = cart_to_spher(rota)

# Calculate the time delays between the Earth center
# and each detector for each sky point
time_delays = [
[
detectors[i].time_delay_from_earth_center(
spher[j][0], spher[j][1], args.trigger_time
)
for j in range(len(spher))
]
for i in range(len(detectors))
]
sky_grid = SkyGrid(
spher[:, 0], spher[:, 1], args.instruments, args.trigger_time
)

extra_attributes = {
'trigger_ra': args.ra,
'trigger_dec': args.dec,
'sky_error': args.sky_error,
'timing_uncertainty': args.timing_uncertainty
}

with HFile(args.output, 'w') as hf:
hf['ra'] = spher[:, 0]
hf['dec'] = spher[:, 1]
hf['trigger_ra'] = [args.ra]
hf['trigger_dec'] = [args.dec]
hf['sky_error'] = [args.sky_error]
hf['trigger_time'] = [args.trigger_time]
hf['timing_uncertainty'] = [args.timing_uncertainty]
hf['instruments'] = [d for d in args.instruments]
hf['time_delays'] = time_delays
sky_grid.write_to_file(args.output, extra_attributes)
12 changes: 11 additions & 1 deletion pycbc/tmpltbank/sky_grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import h5py

from pycbc.detector import Detector
from pycbc.conversions import ensurearray


class SkyGrid:
Expand All @@ -34,6 +35,11 @@ def __init__(self, ra, dec, detectors, ref_gps_time):
# We store the points in a 2D array internally, first dimension runs
# over the list of points, second dimension is RA/dec.
# Question: should we use Astropy sky positions instead?
ra, dec, _ = ensurearray(ra,dec)
if (ra < 0).any() or (ra > 2 * np.pi).any():
raise ValueError('RA must be in the range [0,2π]')
if (dec < -np.pi/2).any() or (dec > np.pi/2).any():
raise ValueError('DEC must be in the range [-π/2, π/2]')
self.positions = np.vstack([ra, dec]).T
self.detectors = sorted(detectors)
self.ref_gps_time = ref_gps_time
Expand Down Expand Up @@ -91,13 +97,17 @@ def read_from_file(cls, path):
ref_gps_time = hf.attrs['ref_gps_time']
return cls(ra, dec, detectors, ref_gps_time)

def write_to_file(self, path):
def write_to_file(self, path, extra_attrs=None, extra_datasets=None):
"""Writes a sky grid to an HDF5 file."""
with h5py.File(path, 'w') as hf:
hf['ra'] = self.ras
hf['dec'] = self.decs
hf.attrs['detectors'] = self.detectors
hf.attrs['ref_gps_time'] = self.ref_gps_time
for attribute in (extra_attrs or {}):
hf.attrs[attribute] = extra_attrs[attribute]
for dataset in (extra_datasets or {}):
hf[dataset] = extra_datasets[dataset]

def calculate_antenna_patterns(self):
"""Calculate the antenna pattern functions at each point in the grid
Expand Down
Loading