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.

Sign up for GitHub

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

Allow RA/dec unit suffix in make_sky_grid #4988

Merged
merged 2 commits into from
Dec 16, 2024
Merged

Allow RA/dec unit suffix in make_sky_grid #4988

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
20b3eeb
Use angle_as_radians in make_skygrid
titodalcanton Dec 13, 2024
47bb4ce
Run black on make_sky_grid
titodalcanton Dec 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
147 changes: 103 additions & 44 deletions bin/pycbc_make_sky_grid
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,15 @@
#!/usr/bin/env python

"""For a given external trigger (GRB, FRB, neutrino, etc...), generate a sky grid covering its localization error region.
"""For a given external trigger (GRB, FRB, neutrino, etc...), generate a sky
grid covering its localization error region.

This sky grid will be used by `pycbc_multi_inspiral` to find multi-detector gravitational wave triggers and calculate the
coherent SNRs and related statistics.
This sky grid will be used by `pycbc_multi_inspiral` to find multi-detector
gravitational wave triggers and calculate the coherent SNRs and related
statistics.

The grid is constructed following the method described in Section V of https://arxiv.org/abs/1410.6042"""
The grid is constructed following the method described in Section V of
https://arxiv.org/abs/1410.6042.
"""

import numpy as np
import argparse
Expand All @@ -15,44 +19,74 @@ 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


def spher_to_cart(sky_points):
"""Convert spherical coordinates to cartesian coordinates.
"""
"""Convert spherical coordinates to cartesian coordinates."""
cart = np.zeros((len(sky_points), 3))
cart[:,0] = np.cos(sky_points[:,0]) * np.cos(sky_points[:,1])
cart[:,1] = np.sin(sky_points[:,0]) * np.cos(sky_points[:,1])
cart[:,2] = np.sin(sky_points[:,1])
cart[:, 0] = np.cos(sky_points[:, 0]) * np.cos(sky_points[:, 1])
cart[:, 1] = np.sin(sky_points[:, 0]) * np.cos(sky_points[:, 1])
cart[:, 2] = np.sin(sky_points[:, 1])
return cart


def cart_to_spher(sky_points):
"""Convert cartesian coordinates to spherical coordinates.
"""
"""Convert cartesian coordinates to spherical coordinates."""
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[:, 0] = np.arctan2(sky_points[:, 1], sky_points[:, 0])
spher[:, 1] = np.arcsin(sky_points[:, 2])
return spher


parser = argparse.ArgumentParser(description=__doc__)
pycbc.add_common_pycbc_options(parser)
parser.add_argument('--ra', type=float,
help="Right ascension (in rad) of the center of the external trigger "
"error box")
parser.add_argument('--dec', type=float,
help="Declination (in rad) of the center of the external trigger "
"error box")
parser.add_argument('--instruments', nargs="+", type=str, required=True,
help="List of instruments to analyze.")
parser.add_argument('--sky-error', type=float, required=True,
help="3-sigma confidence radius (in rad) of the external trigger error "
"box")
parser.add_argument('--trigger-time', type=int, required=True,
help="Time (in s) of the external trigger")
parser.add_argument('--timing-uncertainty', type=float, default=0.0001,
help="Timing uncertainty (in s) we are willing to accept")
parser.add_argument('--output', type=str, required=True,
help="Name of the sky grid")
parser.add_argument(
'--ra',
type=angle_as_radians,
required=True,
help="Right ascension of the center of the external trigger "
"error box. Use the rad or deg suffix to specify units, "
"otherwise radians are assumed.",
)
parser.add_argument(
'--dec',
type=angle_as_radians,
required=True,
help="Declination of the center of the external trigger "
"error box. Use the rad or deg suffix to specify units, "
"otherwise radians are assumed.",
)
parser.add_argument(
'--instruments',
nargs="+",
type=str,
required=True,
help="List of instruments to analyze.",
)
parser.add_argument(
'--sky-error',
type=angle_as_radians,
required=True,
help="3-sigma confidence radius of the external trigger error "
"box. Use the rad or deg suffix to specify units, otherwise "
"radians are assumed.",
)
parser.add_argument(
'--trigger-time',
type=int,
required=True,
help="Time (in s) of the external trigger",
)
parser.add_argument(
'--timing-uncertainty',
type=float,
default=0.0001,
help="Timing uncertainty (in s) we are willing to accept",
)
parser.add_argument(
'--output', type=str, required=True, help="Name of the sky grid"
)

args = parser.parse_args()

Expand All @@ -61,34 +95,51 @@ pycbc.init_logging(args.verbose)
if len(args.instruments) == 1:
parser.error('Can not make a sky grid for only one detector.')

args.instruments.sort() # Put the ifos in alphabetical order
args.instruments.sort() # Put the ifos in alphabetical order
detectors = args.instruments
detectors = [Detector(d) for d in detectors]
detector_pairs = list(itertools.combinations(detectors, 2))

# Calculate the time delay for each detector pair
tds = [detector_pairs[i][0].time_delay_from_detector(detector_pairs[i][1], args.ra, args.dec, args.trigger_time) for i in range(len(detector_pairs))]
tds = [
detector_pairs[i][0].time_delay_from_detector(
detector_pairs[i][1], args.ra, args.dec, args.trigger_time
)
for i in range(len(detector_pairs))
]

# Calculate the light travel time between the detector pairs
light_travel_times = [detector_pairs[i][0].light_travel_time_to_detector(detector_pairs[i][1]) for i in range(len(detector_pairs))]
light_travel_times = [
detector_pairs[i][0].light_travel_time_to_detector(detector_pairs[i][1])
for i in range(len(detector_pairs))
]

# Calculate the required angular spacing between the sky points
ang_spacings = [(2*args.timing_uncertainty) / np.sqrt(light_travel_times[i]**2 - tds[i]**2) for i in range(len(detector_pairs))]
ang_spacings = [
(2 * args.timing_uncertainty)
/ np.sqrt(light_travel_times[i] ** 2 - tds[i] ** 2)
for i in range(len(detector_pairs))
]
angular_spacing = min(ang_spacings)

sky_points = np.zeros((1, 2))

number_of_rings = int(args.sky_error / angular_spacing)

# Generate the sky grid centered at the North pole
for i in range(number_of_rings+1):
for i in range(number_of_rings + 1):
if i == 0:
sky_points[0][0] = 0
sky_points[0][1] = np.pi/2
sky_points[0][1] = np.pi / 2
else:
number_of_points = int(2*np.pi*i)
number_of_points = int(2 * np.pi * i)
for j in range(number_of_points):
sky_points = np.row_stack((sky_points, np.array([j/i, np.pi/2 - i*angular_spacing])))
sky_points = np.row_stack(
(
sky_points,
np.array([j / i, np.pi / 2 - i * angular_spacing]),
)
)

# Convert spherical coordinates to cartesian coordinates
cart = spher_to_cart(sky_points)
Expand All @@ -103,7 +154,7 @@ ort = np.cross(grb_cart, north_pole)
norm = np.linalg.norm(ort)
ort /= norm
n = -np.arccos(np.dot(grb_cart, north_pole))
u = ort*n
u = ort * n

# Rotate the sky grid to the center of the external trigger error box
r = R.from_rotvec(u)
Expand All @@ -112,13 +163,21 @@ 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))]
# 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))
]

with HFile(args.output, 'w') as hf:
hf['ra'] = spher[:,0]
hf['dec'] = spher[:,1]
hf['ra'] = spher[:, 0]
hf['dec'] = spher[:, 1]
hf['trigger_ra'] = [args.ra]
hf['trigger_dec'] = [args.dec]
hf['sky_error'] = [args.sky_error]
Expand Down
Loading