Skip to content

Commit

Permalink
Run black on make_sky_grid
Browse files Browse the repository at this point in the history
  • Loading branch information
titodalcanton committed Dec 13, 2024
1 parent cc7f6b4 commit 6815393
Showing 1 changed file with 99 additions and 44 deletions.
143 changes: 99 additions & 44 deletions bin/pycbc_make_sky_grid
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 @@ -19,44 +23,70 @@ 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=angle_as_radians, required=True,
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,
"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,
"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")
"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 @@ -65,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 @@ -107,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 @@ -116,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

0 comments on commit 6815393

Please sign in to comment.