From 6815393b698bece5a6c713d5378ca8bf87384cae Mon Sep 17 00:00:00 2001 From: Tito Dal Canton Date: Fri, 13 Dec 2024 08:56:59 +0100 Subject: [PATCH] Run black on make_sky_grid --- bin/pycbc_make_sky_grid | 143 +++++++++++++++++++++++++++------------- 1 file changed, 99 insertions(+), 44 deletions(-) diff --git a/bin/pycbc_make_sky_grid b/bin/pycbc_make_sky_grid index 0f3d7fac272..ef2971f9399 100644 --- a/bin/pycbc_make_sky_grid +++ b/bin/pycbc_make_sky_grid @@ -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 @@ -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() @@ -65,19 +95,31 @@ 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)) @@ -85,14 +127,19 @@ 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) @@ -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) @@ -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]