From 99d36dae31e3f66114d2bc3d9a8fec7395ecf4b9 Mon Sep 17 00:00:00 2001 From: Shichao Wu Date: Tue, 31 Oct 2023 16:34:01 +0100 Subject: [PATCH] add coordinates_space.py (#4289) * add coordinates_space.py * add LISA/SSB frame params * add LISA_to_GEO and GEO_to_LISA * add coordinates_space into FieldArray * add doc and Astropy support * update comments on sympy * use fsolve from scipy instead * fix cc issues * fix cc issues * minor fix * update * not use iteration * decouple LISA orbit and more accurate Earth * rename * remove jplephem * add the angular displacement of the Earth * use radians * make func readable in .ini * reverse back to master * correct psi range * reverse to master * fix unit issue in earth_position_SSB * put LISA to the "right" position * add LISA specific transform classes here * change names * update * make a package for coordinates * remove coordinates_space import * move __all__ into __init__.py * remove all coordinates_space * change TIME_OFFSET to seconds * fix SOBHB issue * rename * add SSB or LISA params into fid_params * rename * fix cc issues * fix cc issue * fix cc issue * update * update * fix * add default names * overwrite params with same names * remove pre-fixed names * remove all pre-fixed names * not pop * fix inverse transform * update tc * not overwrite * add SNR support for multi-model * Update waveform.py * t0 issue * t0 issue * Update space.py * add obstime * np.mod(psi_newframe, 2*np.pi) * fix obstime * add support for array inputs * Update hierarchical.py * just use Alex's implementation * CustomTransformMultiOutputs is in another PR, so remove it * add LDC and LAL convention correction * use pycbc standard names * more meaningful name * use pycbc standard names * Update relbin.py * Update parameters.py * remove unnecessary changes * fix cc issue * fix cc issue * fix cc issue * fix cc issue * compactify * compactify * add __all__ back * Update transforms.py * Update transforms.py * Update test_transforms.py * Update transforms.py * update doc * fix time warning * Update space.py * Update test_transforms.py * Create test_coordinates_space.py * fix cc issues * fix cc issues * fix cc issue * Update tox.ini * Update tox.ini * Update tox.ini * Update tox.ini * Update tox.ini * Update tox.ini * Update tox.ini * Update test_coordinates_space.py * add inline doc * Update tox.ini * add check of bbhx * Update test_coordinates_space.py * Update tox.ini * Update test_coordinates_space.py * add MultibandRelativeTimeDom into hierarchical.py * Update __init__.py * Update hierarchical.py * Update hierarchical.py * Update relbin.py * Update hierarchical.py * Update hierarchical.py * Update relbin.py * Update hierarchical.py * Update hierarchical.py * Update hierarchical.py * Update hierarchical.py * Update hierarchical.py * Update relbin.py * Update hierarchical.py * Update hierarchical.py * Update relbin.py * Update __init__.py * Update space.py * Update space.py * Update space.py * fix psi issue * Update test_coordinates_space.py * Update test_coordinates_space.py --- pycbc/coordinates/__init__.py | 37 + pycbc/{coordinates.py => coordinates/base.py} | 17 +- pycbc/coordinates/space.py | 949 ++++++++++++++++++ pycbc/transforms.py | 519 ++++++++++ test/test_coordinates_space.py | 358 +++++++ test/test_transforms.py | 7 +- tox.ini | 24 +- 7 files changed, 1905 insertions(+), 6 deletions(-) create mode 100644 pycbc/coordinates/__init__.py rename pycbc/{coordinates.py => coordinates/base.py} (92%) create mode 100644 pycbc/coordinates/space.py create mode 100644 test/test_coordinates_space.py diff --git a/pycbc/coordinates/__init__.py b/pycbc/coordinates/__init__.py new file mode 100644 index 00000000000..e1d7c81d90d --- /dev/null +++ b/pycbc/coordinates/__init__.py @@ -0,0 +1,37 @@ +# Copyright (C) 2023 Shichao Wu, Alex Nitz +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +""" +This modules provides functions for base coordinate transformations, and +more advanced transformations between ground-based detectors and space-borne +detectors. +""" + +from pycbc.coordinates.base import * +from pycbc.coordinates.space import * + + +__all__ = ['cartesian_to_spherical_rho', 'cartesian_to_spherical_azimuthal', + 'cartesian_to_spherical_polar', 'cartesian_to_spherical', + 'spherical_to_cartesian', + 'TIME_OFFSET_20_DEGREES', + 'localization_to_propagation_vector', + 'propagation_vector_to_localization', 'polarization_newframe', + 't_lisa_from_ssb', 't_ssb_from_t_lisa', + 'ssb_to_lisa', 'lisa_to_ssb', + 'rotation_matrix_ssb_to_lisa', 'rotation_matrix_ssb_to_geo', + 'lisa_position_ssb', 'earth_position_ssb', + 't_geo_from_ssb', 't_ssb_from_t_geo', 'ssb_to_geo', 'geo_to_ssb', + 'lisa_to_geo', 'geo_to_lisa', + ] diff --git a/pycbc/coordinates.py b/pycbc/coordinates/base.py similarity index 92% rename from pycbc/coordinates.py rename to pycbc/coordinates/base.py index dab734b1aa5..48dffe8573b 100644 --- a/pycbc/coordinates.py +++ b/pycbc/coordinates/base.py @@ -14,7 +14,18 @@ # You should have received a copy of the GNU General Public License along # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -""" Coordinate transformations. + + +# +# ============================================================================= +# +# Preamble +# +# ============================================================================= +# +""" +Base coordinate transformations, this module provides transformations between +cartesian and spherical coordinates. """ import numpy @@ -60,6 +71,7 @@ def cartesian_to_spherical_azimuthal(x, y): phi = numpy.arctan2(y, x) return phi % (2 * numpy.pi) + def cartesian_to_spherical_polar(x, y, z): """ Calculates the polar angle in spherical coordinates from Cartesian coordinates. The polar angle is in [0,pi]. @@ -141,7 +153,8 @@ def spherical_to_cartesian(rho, phi, theta): z = rho * numpy.cos(theta) return x, y, z + __all__ = ['cartesian_to_spherical_rho', 'cartesian_to_spherical_azimuthal', 'cartesian_to_spherical_polar', 'cartesian_to_spherical', 'spherical_to_cartesian', - ] + ] diff --git a/pycbc/coordinates/space.py b/pycbc/coordinates/space.py new file mode 100644 index 00000000000..4a027746b0b --- /dev/null +++ b/pycbc/coordinates/space.py @@ -0,0 +1,949 @@ +# Copyright (C) 2023 Shichao Wu, Alex Nitz +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + + +# +# ============================================================================= +# +# Preamble +# +# ============================================================================= +# +""" +This module provides coordinate transformations related to space-borne +detectors, such as coordinate transformations between space-borne detectors +and ground-based detectors. Note that current LISA orbit used in this module +is a circular orbit, need to be replaced by a more realistic and general orbit +model in the near future. +""" + +import numpy as np +from scipy.spatial.transform import Rotation +from scipy.optimize import fsolve +from astropy import units +from astropy.constants import c, au +from astropy.time import Time +from astropy.coordinates import BarycentricMeanEcliptic, PrecessedGeocentric +from astropy.coordinates import get_body_barycentric +from astropy.coordinates import SkyCoord +from astropy.coordinates.builtin_frames import ecliptic_transforms + +# This constant makes sure LISA is behind the Earth by 19-23 degrees. +# Making this a stand-alone constant will also make it callable by +# the waveform plugin and PE config file. In the unit of 's'. +TIME_OFFSET_20_DEGREES = 7365189.431698299 + +# "rotation_matrix_ssb_to_lisa" and "lisa_position_ssb" should be +# more general for other detectors in the near future. + + +def rotation_matrix_ssb_to_lisa(alpha): + """ The rotation matrix (of frame basis) from SSB frame to LISA frame. + This function assumes the angle between LISA plane and the ecliptic + is 60 degrees, and the period of LISA's self-rotation and orbital + revolution is both one year. + + Parameters + ---------- + alpha : float + The angular displacement of LISA in SSB frame. + In the unit of 'radian'. + + Returns + ------- + r_total : numpy.array + A 3x3 rotation matrix from SSB frame to LISA frame. + """ + r = Rotation.from_rotvec([ + [0, 0, alpha], + [0, -np.pi/3, 0], + [0, 0, -alpha] + ]).as_matrix() + r_total = np.array(r[0]) @ np.array(r[1]) @ np.array(r[2]) + + return r_total + + +def lisa_position_ssb(t_lisa, t0=TIME_OFFSET_20_DEGREES): + """ Calculating the position vector and angular displacement of LISA + in the SSB frame, at a given time. This function assumes LISA's barycenter + is orbiting around a circular orbit within the ecliptic behind the Earth. + The period of it is one year. + + Parameters + ---------- + t_lisa : float + The time when a GW signal arrives at the origin of LISA frame, + or any other time you want. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + + Returns + ------- + (p, alpha) : tuple + p : numpy.array + The position vector of LISA in the SSB frame. In the unit of 'm'. + alpha : float + The angular displacement of LISA in the SSB frame. + In the unit of 'radian'. + """ + OMEGA_0 = 1.99098659277e-7 + R_ORBIT = au.value + alpha = np.mod(OMEGA_0 * (t_lisa + t0), 2*np.pi) + p = np.array([[R_ORBIT * np.cos(alpha)], + [R_ORBIT * np.sin(alpha)], + [0]], dtype=object) + return (p, alpha) + + +def localization_to_propagation_vector(longitude, latitude, + use_astropy=True, frame=None): + """ Converting the sky localization to the corresponding + propagation unit vector of a GW signal. + + Parameters + ---------- + longitude : float + The longitude, in the unit of 'radian'. + latitude : float + The latitude, in the unit of 'radian'. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + frame : astropy.coordinates + The frame from astropy.coordinates if use_astropy is True, + the default is None. + + Returns + ------- + [[x], [y], [z]] : numpy.array + The propagation unit vector of that GW signal. + """ + if use_astropy: + x = -frame.cartesian.x.value + y = -frame.cartesian.y.value + z = -frame.cartesian.z.value + else: + x = -np.cos(latitude) * np.cos(longitude) + y = -np.cos(latitude) * np.sin(longitude) + z = -np.sin(latitude) + v = np.array([[x], [y], [z]]) + + return v / np.linalg.norm(v) + + +def propagation_vector_to_localization(k, use_astropy=True, frame=None): + """ Converting the propagation unit vector to the corresponding + sky localization of a GW signal. + + Parameters + ---------- + k : numpy.array + The propagation unit vector of a GW signal. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + frame : astropy.coordinates + The frame from astropy.coordinates if use_astropy is True, + the default is None. + + Returns + ------- + (longitude, latitude) : tuple + The sky localization of that GW signal. + """ + if use_astropy: + try: + longitude = frame.lon.rad + latitude = frame.lat.rad + except AttributeError: + longitude = frame.ra.rad + latitude = frame.dec.rad + else: + # latitude already within [-pi/2, pi/2] + latitude = np.float64(np.arcsin(-k[2])) + longitude = np.float64(np.arctan2(-k[1]/np.cos(latitude), + -k[0]/np.cos(latitude))) + # longitude should within [0, 2*pi) + longitude = np.mod(longitude, 2*np.pi) + + return (longitude, latitude) + + +def polarization_newframe(polarization, k, rotation_matrix, use_astropy=True, + old_frame=None, new_frame=None): + """ Converting a polarization angle from a frame to a new frame + by using rotation matrix method. + + Parameters + ---------- + polarization : float + The polarization angle in the old frame, in the unit of 'radian'. + k : numpy.array + The propagation unit vector of a GW signal in the old frame. + rotation_matrix : numpy.array + The rotation matrix (of frame basis) from the old frame to + the new frame. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + old_frame : astropy.coordinates + The frame from astropy.coordinates if use_astropy is True, + the default is None. + new_frame : astropy.coordinates + The frame from astropy.coordinates if use_astropy is True, + the default is None. The new frame for the new polarization + angle. + + Returns + ------- + polarization_new_frame : float + The polarization angle in the new frame of that GW signal. + """ + longitude, _ = propagation_vector_to_localization( + k, use_astropy, old_frame) + u = np.array([[np.sin(longitude)], [-np.cos(longitude)], [0]]) + rotation_vector = polarization * k + rotation_polarization = Rotation.from_rotvec(rotation_vector.T[0]) + p = rotation_polarization.apply(u.T[0]).reshape(3, 1) + p_newframe = rotation_matrix.T @ p + k_newframe = rotation_matrix.T @ k + longitude_newframe, latitude_newframe = \ + propagation_vector_to_localization(k_newframe, use_astropy, new_frame) + u_newframe = np.array([[np.sin(longitude_newframe)], + [-np.cos(longitude_newframe)], [0]]) + v_newframe = np.array([ + [-np.sin(latitude_newframe) * np.cos(longitude_newframe)], + [-np.sin(latitude_newframe) * np.sin(longitude_newframe)], + [np.cos(latitude_newframe)]]) + p_dot_u_newframe = np.vdot(p_newframe, u_newframe) + p_dot_v_newframe = np.vdot(p_newframe, v_newframe) + polarization_new_frame = np.arctan2(p_dot_v_newframe, p_dot_u_newframe) + polarization_new_frame = np.mod(polarization_new_frame, 2*np.pi) + # avoid the round error + if polarization_new_frame == 2*np.pi: + polarization_new_frame = 0 + + return polarization_new_frame + + +def t_lisa_from_ssb(t_ssb, longitude_ssb, latitude_ssb, + t0=TIME_OFFSET_20_DEGREES): + """ Calculating the time when a GW signal arrives at the barycenter + of LISA, by using the time and sky localization in SSB frame. + + Parameters + ---------- + t_ssb : float + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + + Returns + ------- + t_lisa : float + The time when a GW signal arrives at the origin of LISA frame. + """ + k = localization_to_propagation_vector( + longitude_ssb, latitude_ssb, use_astropy=False) + + def equation(t_lisa): + # LISA is moving, when GW arrives at LISA center, + # time is t_lisa, not t_ssb. + p = lisa_position_ssb(t_lisa, t0)[0] + return t_lisa - t_ssb - np.vdot(k, p) / c.value + + return fsolve(equation, t_ssb)[0] + + +def t_ssb_from_t_lisa(t_lisa, longitude_ssb, latitude_ssb, + t0=TIME_OFFSET_20_DEGREES): + """ Calculating the time when a GW signal arrives at the barycenter + of SSB, by using the time in LISA frame and sky localization in SSB frame. + + Parameters + ---------- + t_lisa : float + The time when a GW signal arrives at the origin of LISA frame. + In the unit of 's'. + longitude_ssb : float + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + + Returns + ------- + t_ssb : float + The time when a GW signal arrives at the origin of SSB frame. + """ + k = localization_to_propagation_vector( + longitude_ssb, latitude_ssb, use_astropy=False) + # LISA is moving, when GW arrives at LISA center, + # time is t_lisa, not t_ssb. + p = lisa_position_ssb(t_lisa, t0)[0] + + def equation(t_ssb): + return t_lisa - t_ssb - np.vdot(k, p) / c.value + + return fsolve(equation, t_lisa)[0] + + +def ssb_to_lisa(t_ssb, longitude_ssb, latitude_ssb, polarization_ssb, + t0=TIME_OFFSET_20_DEGREES): + """ Converting the arrive time, the sky localization, and the polarization + from the SSB frame to the LISA frame. + + Parameters + ---------- + t_ssb : float or numpy.array + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float or numpy.array + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float or numpy.array + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + polarization_ssb : float or numpy.array + The polarization angle of a GW signal in SSB frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + + Returns + ------- + (t_lisa, longitude_lisa, latitude_lisa, polarization_lisa) : tuple + t_lisa : float or numpy.array + The time when a GW signal arrives at the origin of LISA frame. + In the unit of 's'. + longitude_lisa : float or numpy.array + The longitude of a GW signal in LISA frame, in the unit of 'radian'. + latitude_lisa : float or numpy.array + The latitude of a GW signal in LISA frame, in the unit of 'radian'. + polarization_lisa : float or numpy.array + The polarization angle of a GW signal in LISA frame. + In the unit of 'radian'. + """ + if not isinstance(t_ssb, np.ndarray): + t_ssb = np.array([t_ssb]) + if not isinstance(longitude_ssb, np.ndarray): + longitude_ssb = np.array([longitude_ssb]) + if not isinstance(latitude_ssb, np.ndarray): + latitude_ssb = np.array([latitude_ssb]) + if not isinstance(polarization_ssb, np.ndarray): + polarization_ssb = np.array([polarization_ssb]) + num = len(t_ssb) + t_lisa, longitude_lisa = np.zeros(num), np.zeros(num) + latitude_lisa, polarization_lisa = np.zeros(num), np.zeros(num) + + for i in range(num): + if longitude_ssb[i] < 0 or longitude_ssb[i] >= 2*np.pi: + raise ValueError("Longitude should within [0, 2*pi).") + if latitude_ssb[i] < -np.pi/2 or latitude_ssb[i] > np.pi/2: + raise ValueError("Latitude should within [-pi/2, pi/2].") + if polarization_ssb[i] < 0 or polarization_ssb[i] >= 2*np.pi: + raise ValueError("Polarization angle should within [0, 2*pi).") + t_lisa[i] = t_lisa_from_ssb(t_ssb[i], longitude_ssb[i], + latitude_ssb[i], t0) + k_ssb = localization_to_propagation_vector( + longitude_ssb[i], latitude_ssb[i], use_astropy=False) + # Although t_lisa calculated above using the corrected LISA position + # vector by adding t0, it corresponds to the true t_ssb, not t_ssb+t0, + # we need to include t0 again to correct LISA position. + alpha = lisa_position_ssb(t_lisa[i], t0)[1] + rotation_matrix_lisa = rotation_matrix_ssb_to_lisa(alpha) + k_lisa = rotation_matrix_lisa.T @ k_ssb + longitude_lisa[i], latitude_lisa[i] = \ + propagation_vector_to_localization(k_lisa, use_astropy=False) + polarization_lisa[i] = polarization_newframe( + polarization_ssb[i], k_ssb, rotation_matrix_lisa, + use_astropy=False) + + if num == 1: + params_lisa = (t_lisa[0], longitude_lisa[0], + latitude_lisa[0], polarization_lisa[0]) + else: + params_lisa = (t_lisa, longitude_lisa, + latitude_lisa, polarization_lisa) + + return params_lisa + + +def lisa_to_ssb(t_lisa, longitude_lisa, latitude_lisa, polarization_lisa, + t0=TIME_OFFSET_20_DEGREES): + """ Converting the arrive time, the sky localization, and the polarization + from the LISA frame to the SSB frame. + + Parameters + ---------- + t_lisa : float or numpy.array + The time when a GW signal arrives at the origin of LISA frame. + In the unit of 's'. + longitude_lisa : float or numpy.array + The longitude of a GW signal in LISA frame, in the unit of 'radian'. + latitude_lisa : float or numpy.array + The latitude of a GW signal in LISA frame, in the unit of 'radian'. + polarization_lisa : float or numpy.array + The polarization angle of a GW signal in LISA frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + + Returns + ------- + (t_ssb, longitude_ssb, latitude_ssb, polarization_ssb) : tuple + t_ssb : float or numpy.array + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float or numpy.array + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float or numpy.array + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + polarization_ssb : float or numpy.array + The polarization angle of a GW signal in SSB frame. + In the unit of 'radian'. + """ + if not isinstance(t_lisa, np.ndarray): + t_lisa = np.array([t_lisa]) + if not isinstance(longitude_lisa, np.ndarray): + longitude_lisa = np.array([longitude_lisa]) + if not isinstance(latitude_lisa, np.ndarray): + latitude_lisa = np.array([latitude_lisa]) + if not isinstance(polarization_lisa, np.ndarray): + polarization_lisa = np.array([polarization_lisa]) + num = len(t_lisa) + t_ssb, longitude_ssb = np.zeros(num), np.zeros(num) + latitude_ssb, polarization_ssb = np.zeros(num), np.zeros(num) + + for i in range(num): + if longitude_lisa[i] < 0 or longitude_lisa[i] >= 2*np.pi: + raise ValueError("Longitude should within [0, 2*pi).") + if latitude_lisa[i] < -np.pi/2 or latitude_lisa[i] > np.pi/2: + raise ValueError("Latitude should within [-pi/2, pi/2].") + if polarization_lisa[i] < 0 or polarization_lisa[i] >= 2*np.pi: + raise ValueError("Polarization angle should within [0, 2*pi).") + k_lisa = localization_to_propagation_vector( + longitude_lisa[i], latitude_lisa[i], use_astropy=False) + alpha = lisa_position_ssb(t_lisa[i], t0)[1] + rotation_matrix_lisa = rotation_matrix_ssb_to_lisa(alpha) + k_ssb = rotation_matrix_lisa @ k_lisa + longitude_ssb[i], latitude_ssb[i] = \ + propagation_vector_to_localization(k_ssb, use_astropy=False) + t_ssb[i] = t_ssb_from_t_lisa(t_lisa[i], longitude_ssb[i], + latitude_ssb[i], t0) + polarization_ssb[i] = polarization_newframe( + polarization_lisa[i], k_lisa, rotation_matrix_lisa.T, + use_astropy=False) + + if num == 1: + params_ssb = (t_ssb[0], longitude_ssb[0], + latitude_ssb[0], polarization_ssb[0]) + else: + params_ssb = (t_ssb, longitude_ssb, + latitude_ssb, polarization_ssb) + + return params_ssb + + +def rotation_matrix_ssb_to_geo(epsilon=np.deg2rad(23.439281)): + """ The rotation matrix (of frame basis) from SSB frame to + geocentric frame. + + Parameters + ---------- + epsilon : float + The Earth's axial tilt (obliquity), in the unit of 'radian'. + + Returns + ------- + r : numpy.array + A 3x3 rotation matrix from SSB frame to geocentric frame. + """ + r = Rotation.from_rotvec([ + [-epsilon, 0, 0] + ]).as_matrix() + + return np.array(r[0]) + + +def earth_position_ssb(t_geo): + """ Calculating the position vector and angular displacement of the Earth + in the SSB frame, at a given time. By using Astropy. + + Parameters + ---------- + t_geo : float + The time when a GW signal arrives at the origin of geocentric frame, + or any other time you want. + + Returns + ------- + (p, alpha) : tuple + p : numpy.array + The position vector of the Earth in the SSB frame. In the unit of 'm'. + alpha : float + The angular displacement of the Earth in the SSB frame. + In the unit of 'radian'. + """ + t = Time(t_geo, format='gps') + pos = get_body_barycentric('earth', t) + # BarycentricMeanEcliptic doesn't have obstime attribute, + # it's a good inertial frame, but ICRS is not. + icrs_coord = SkyCoord(pos, frame='icrs', obstime=t) + bme_coord = icrs_coord.transform_to( + BarycentricMeanEcliptic(equinox='J2000')) + x = bme_coord.cartesian.x.to(units.m).value + y = bme_coord.cartesian.y.to(units.m).value + z = bme_coord.cartesian.z.to(units.m).value + p = np.array([[x], [y], [z]]) + alpha = bme_coord.lon.rad + + return (p, alpha) + + +def t_geo_from_ssb(t_ssb, longitude_ssb, latitude_ssb, + use_astropy=True, frame=None): + """ Calculating the time when a GW signal arrives at the barycenter + of the Earth, by using the time and sky localization in SSB frame. + + Parameters + ---------- + t_ssb : float + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + + Returns + ------- + t_geo : float + The time when a GW signal arrives at the origin of geocentric frame. + """ + k = localization_to_propagation_vector( + longitude_ssb, latitude_ssb, use_astropy, frame) + + def equation(t_geo): + # Earth is moving, when GW arrives at Earth center, + # time is t_geo, not t_ssb. + p = earth_position_ssb(t_geo)[0] + return t_geo - t_ssb - np.vdot(k, p) / c.value + + return fsolve(equation, t_ssb)[0] + + +def t_ssb_from_t_geo(t_geo, longitude_ssb, latitude_ssb, + use_astropy=True, frame=None): + """ Calculating the time when a GW signal arrives at the barycenter + of SSB, by using the time in geocentric frame and sky localization + in SSB frame. + + Parameters + ---------- + t_geo : float + The time when a GW signal arrives at the origin of geocentric frame. + In the unit of 's'. + longitude_ssb : float + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + + Returns + ------- + t_ssb : float + The time when a GW signal arrives at the origin of SSB frame. + """ + k = localization_to_propagation_vector( + longitude_ssb, latitude_ssb, use_astropy, frame) + # Earth is moving, when GW arrives at Earth center, + # time is t_geo, not t_ssb. + p = earth_position_ssb(t_geo)[0] + + def equation(t_ssb): + return t_geo - t_ssb - np.vdot(k, p) / c.value + + return fsolve(equation, t_geo)[0] + + +def ssb_to_geo(t_ssb, longitude_ssb, latitude_ssb, polarization_ssb, + use_astropy=True): + """ Converting the arrive time, the sky localization, and the polarization + from the SSB frame to the geocentric frame. + + Parameters + ---------- + t_ssb : float or numpy.array + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float or numpy.array + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float or numpy.array + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + polarization_ssb : float or numpy.array + The polarization angle of a GW signal in SSB frame. + In the unit of 'radian'. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + + Returns + ------- + (t_geo, longitude_geo, latitude_geo, polarization_geo) : tuple + t_geo : float or numpy.array + The time when a GW signal arrives at the origin of geocentric frame. + In the unit of 's'. + longitude_geo : float or numpy.array + The longitude of a GW signal in geocentric frame. + In the unit of 'radian'. + latitude_geo : float or numpy.array + The latitude of a GW signal in geocentric frame. + In the unit of 'radian'. + polarization_geo : float or numpy.array + The polarization angle of a GW signal in geocentric frame. + In the unit of 'radian'. + """ + if not isinstance(t_ssb, np.ndarray): + t_ssb = np.array([t_ssb]) + if not isinstance(longitude_ssb, np.ndarray): + longitude_ssb = np.array([longitude_ssb]) + if not isinstance(latitude_ssb, np.ndarray): + latitude_ssb = np.array([latitude_ssb]) + if not isinstance(polarization_ssb, np.ndarray): + polarization_ssb = np.array([polarization_ssb]) + num = len(t_ssb) + t_geo = np.full(num, np.nan) + longitude_geo = np.full(num, np.nan) + latitude_geo = np.full(num, np.nan) + polarization_geo = np.full(num, np.nan) + + for i in range(num): + if longitude_ssb[i] < 0 or longitude_ssb[i] >= 2*np.pi: + raise ValueError("Longitude should within [0, 2*pi).") + if latitude_ssb[i] < -np.pi/2 or latitude_ssb[i] > np.pi/2: + raise ValueError("Latitude should within [-pi/2, pi/2].") + if polarization_ssb[i] < 0 or polarization_ssb[i] >= 2*np.pi: + raise ValueError("Polarization angle should within [0, 2*pi).") + + if use_astropy: + # BarycentricMeanEcliptic doesn't have obstime attribute, + # it's a good inertial frame, but PrecessedGeocentric is not. + bme_coord = BarycentricMeanEcliptic( + lon=longitude_ssb[i]*units.radian, + lat=latitude_ssb[i]*units.radian, + equinox='J2000') + t_geo[i] = t_geo_from_ssb(t_ssb[i], longitude_ssb[i], + latitude_ssb[i], use_astropy, bme_coord) + geo_sky = bme_coord.transform_to(PrecessedGeocentric( + equinox='J2000', obstime=Time(t_geo[i], format='gps'))) + longitude_geo[i] = geo_sky.ra.rad + latitude_geo[i] = geo_sky.dec.rad + k_geo = localization_to_propagation_vector( + longitude_geo[i], latitude_geo[i], + use_astropy, geo_sky) + k_ssb = localization_to_propagation_vector( + None, None, use_astropy, bme_coord) + rotation_matrix_geo = \ + ecliptic_transforms.icrs_to_baryecliptic( + from_coo=None, + to_frame=BarycentricMeanEcliptic(equinox='J2000')) + polarization_geo[i] = polarization_newframe( + polarization_ssb[i], k_ssb, + rotation_matrix_geo, use_astropy, + old_frame=bme_coord, + new_frame=geo_sky) + else: + t_geo[i] = t_geo_from_ssb(t_ssb[i], longitude_ssb[i], + latitude_ssb[i], use_astropy) + rotation_matrix_geo = rotation_matrix_ssb_to_geo() + k_ssb = localization_to_propagation_vector( + longitude_ssb[i], latitude_ssb[i], + use_astropy) + k_geo = rotation_matrix_geo.T @ k_ssb + longitude_geo[i], latitude_geo[i] = \ + propagation_vector_to_localization(k_geo, use_astropy) + polarization_geo[i] = polarization_newframe( + polarization_ssb[i], k_ssb, + rotation_matrix_geo, use_astropy) + + # As mentioned in LDC manual, the p,q vectors are opposite between + # LDC and LAL conventions, see Sec 4.1.5 in . + polarization_geo[i] = np.mod(polarization_geo[i]+np.pi, 2*np.pi) + + if num == 1: + params_geo = (t_geo[0], longitude_geo[0], + latitude_geo[0], polarization_geo[0]) + else: + params_geo = (t_geo, longitude_geo, + latitude_geo, polarization_geo) + + return params_geo + + +def geo_to_ssb(t_geo, longitude_geo, latitude_geo, polarization_geo, + use_astropy=True): + """ Converting the arrive time, the sky localization, and the polarization + from the geocentric frame to the SSB frame. + + Parameters + ---------- + t_geo : float or numpy.array + The time when a GW signal arrives at the origin of geocentric frame. + In the unit of 's'. + longitude_geo : float or numpy.array + The longitude of a GW signal in geocentric frame. + In the unit of 'radian'. + latitude_geo : float or numpy.array + The latitude of a GW signal in geocentric frame. + In the unit of 'radian'. + polarization_geo : float or numpy.array + The polarization angle of a GW signal in geocentric frame. + In the unit of 'radian'. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + + Returns + ------- + (t_ssb, longitude_ssb, latitude_ssb, polarization_ssb) : tuple + t_ssb : float or numpy.array + The time when a GW signal arrives at the origin of SSB frame. + In the unit of 's'. + longitude_ssb : float or numpy.array + The ecliptic longitude of a GW signal in SSB frame. + In the unit of 'radian'. + latitude_ssb : float or numpy.array + The ecliptic latitude of a GW signal in SSB frame. + In the unit of 'radian'. + polarization_ssb : float or numpy.array + The polarization angle of a GW signal in SSB frame. + In the unit of 'radian'. + """ + if not isinstance(t_geo, np.ndarray): + t_geo = np.array([t_geo]) + if not isinstance(longitude_geo, np.ndarray): + longitude_geo = np.array([longitude_geo]) + if not isinstance(latitude_geo, np.ndarray): + latitude_geo = np.array([latitude_geo]) + if not isinstance(polarization_geo, np.ndarray): + polarization_geo = np.array([polarization_geo]) + num = len(t_geo) + t_ssb = np.full(num, np.nan) + longitude_ssb = np.full(num, np.nan) + latitude_ssb = np.full(num, np.nan) + polarization_ssb = np.full(num, np.nan) + + for i in range(num): + if longitude_geo[i] < 0 or longitude_geo[i] >= 2*np.pi: + raise ValueError("Longitude should within [0, 2*pi).") + if latitude_geo[i] < -np.pi/2 or latitude_geo[i] > np.pi/2: + raise ValueError("Latitude should within [-pi/2, pi/2].") + if polarization_geo[i] < 0 or polarization_geo[i] >= 2*np.pi: + raise ValueError("Polarization angle should within [0, 2*pi).") + + if use_astropy: + # BarycentricMeanEcliptic doesn't have obstime attribute, + # it's a good inertial frame, but PrecessedGeocentric is not. + geo_coord = PrecessedGeocentric( + ra=longitude_geo[i]*units.radian, + dec=latitude_geo[i]*units.radian, + equinox='J2000', + obstime=Time(t_geo[i], format='gps')) + ssb_sky = geo_coord.transform_to( + BarycentricMeanEcliptic(equinox='J2000')) + longitude_ssb[i] = ssb_sky.lon.rad + latitude_ssb[i] = ssb_sky.lat.rad + k_ssb = localization_to_propagation_vector( + longitude_ssb[i], latitude_ssb[i], + use_astropy, ssb_sky) + k_geo = localization_to_propagation_vector( + None, None, use_astropy, geo_coord) + rotation_matrix_geo = \ + ecliptic_transforms.icrs_to_baryecliptic( + from_coo=None, + to_frame=BarycentricMeanEcliptic(equinox='J2000')) + t_ssb[i] = t_ssb_from_t_geo(t_geo[i], longitude_ssb[i], + latitude_ssb[i], use_astropy, + ssb_sky) + polarization_ssb[i] = polarization_newframe( + polarization_geo[i], k_geo, + rotation_matrix_geo.T, + use_astropy, + old_frame=geo_coord, + new_frame=ssb_sky) + else: + rotation_matrix_geo = rotation_matrix_ssb_to_geo() + k_geo = localization_to_propagation_vector( + longitude_geo[i], latitude_geo[i], use_astropy) + k_ssb = rotation_matrix_geo @ k_geo + longitude_ssb[i], latitude_ssb[i] = \ + propagation_vector_to_localization(k_ssb, use_astropy) + t_ssb[i] = t_ssb_from_t_geo(t_geo[i], longitude_ssb[i], + latitude_ssb[i], use_astropy) + polarization_ssb[i] = polarization_newframe( + polarization_geo[i], k_geo, + rotation_matrix_geo.T, use_astropy) + + # As mentioned in LDC manual, the p,q vectors are opposite between + # LDC and LAL conventions, see Sec 4.1.5 in . + polarization_ssb[i] = np.mod(polarization_ssb[i]-np.pi, 2*np.pi) + + if num == 1: + params_ssb = (t_ssb[0], longitude_ssb[0], + latitude_ssb[0], polarization_ssb[0]) + else: + params_ssb = (t_ssb, longitude_ssb, + latitude_ssb, polarization_ssb) + + return params_ssb + + +def lisa_to_geo(t_lisa, longitude_lisa, latitude_lisa, polarization_lisa, + t0=TIME_OFFSET_20_DEGREES, use_astropy=True): + """ Converting the arrive time, the sky localization, and the polarization + from the LISA frame to the geocentric frame. + + Parameters + ---------- + t_lisa : float or numpy.array + The time when a GW signal arrives at the origin of LISA frame. + In the unit of 's'. + longitude_lisa : float or numpy.array + The longitude of a GW signal in LISA frame, in the unit of 'radian'. + latitude_lisa : float or numpy.array + The latitude of a GW signal in LISA frame, in the unit of 'radian'. + polarization_lisa : float or numpy.array + The polarization angle of a GW signal in LISA frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + + Returns + ------- + (t_geo, longitude_geo, latitude_geo, polarization_geo) : tuple + t_geo : float or numpy.array + The time when a GW signal arrives at the origin of geocentric frame. + In the unit of 's'. + longitude_geo : float or numpy.array + The ecliptic longitude of a GW signal in geocentric frame. + In the unit of 'radian'. + latitude_geo : float or numpy.array + The ecliptic latitude of a GW signal in geocentric frame. + In the unit of 'radian'. + polarization_geo : float or numpy.array + The polarization angle of a GW signal in geocentric frame. + In the unit of 'radian'. + """ + t_ssb, longitude_ssb, latitude_ssb, polarization_ssb = lisa_to_ssb( + t_lisa, longitude_lisa, latitude_lisa, polarization_lisa, t0) + t_geo, longitude_geo, latitude_geo, polarization_geo = ssb_to_geo( + t_ssb, longitude_ssb, latitude_ssb, polarization_ssb, use_astropy) + + return (t_geo, longitude_geo, latitude_geo, polarization_geo) + + +def geo_to_lisa(t_geo, longitude_geo, latitude_geo, polarization_geo, + t0=TIME_OFFSET_20_DEGREES, use_astropy=True): + """ Converting the arrive time, the sky localization, and the polarization + from the geocentric frame to the LISA frame. + + Parameters + ---------- + t_geo : float or numpy.array + The time when a GW signal arrives at the origin of geocentric frame. + In the unit of 's'. + longitude_geo : float or numpy.array + The longitude of a GW signal in geocentric frame. + In the unit of 'radian'. + latitude_geo : float or numpy.array + The latitude of a GW signal in geocentric frame. + In the unit of 'radian'. + polarization_geo : float or numpy.array + The polarization angle of a GW signal in geocentric frame. + In the unit of 'radian'. + t0 : float + The initial time offset of LISA, in the unit of 's', + default is 7365189.431698299. This makes sure LISA is behind + the Earth by 19-23 degrees. + use_astropy : bool + Using Astropy to calculate the sky localization or not. + Default is True. + + Returns + ------- + (t_lisa, longitude_lisa, latitude_lisa, polarization_lisa) : tuple + t_lisa : float or numpy.array + The time when a GW signal arrives at the origin of LISA frame. + In the unit of 's'. + longitude_lisa : float or numpy.array + The longitude of a GW signal in LISA frame, in the unit of 'radian'. + latitude_lisa : float or numpy.array + The latitude of a GW signal in LISA frame, in the unit of 'radian'. + polarization_geo : float or numpy.array + The polarization angle of a GW signal in LISA frame. + In the unit of 'radian'. + """ + t_ssb, longitude_ssb, latitude_ssb, polarization_ssb = geo_to_ssb( + t_geo, longitude_geo, latitude_geo, polarization_geo, use_astropy) + t_lisa, longitude_lisa, latitude_lisa, polarization_lisa = ssb_to_lisa( + t_ssb, longitude_ssb, latitude_ssb, polarization_ssb, t0) + + return (t_lisa, longitude_lisa, latitude_lisa, polarization_lisa) + + +__all__ = ['TIME_OFFSET_20_DEGREES', + 'localization_to_propagation_vector', + 'propagation_vector_to_localization', 'polarization_newframe', + 't_lisa_from_ssb', 't_ssb_from_t_lisa', + 'ssb_to_lisa', 'lisa_to_ssb', + 'rotation_matrix_ssb_to_lisa', 'rotation_matrix_ssb_to_geo', + 'lisa_position_ssb', 'earth_position_ssb', + 't_geo_from_ssb', 't_ssb_from_t_geo', 'ssb_to_geo', 'geo_to_ssb', + 'lisa_to_geo', 'geo_to_lisa', + ] diff --git a/pycbc/transforms.py b/pycbc/transforms.py index d4e9f79894e..c8fa44a9a9b 100644 --- a/pycbc/transforms.py +++ b/pycbc/transforms.py @@ -1468,6 +1468,402 @@ def from_config(cls, cp, section, outputs): ) +class GEOToSSB(BaseTransform): + """Converts arrival time, sky localization, and polarization angle in the + geocentric frame to the corresponding values in the SSB frame.""" + + name = "geo_to_ssb" + + default_params_name = { + 'default_tc_geo': parameters.tc, + 'default_longitude_geo': parameters.ra, + 'default_latitude_geo': parameters.dec, + 'default_polarization_geo': parameters.polarization, + 'default_tc_ssb': parameters.tc, + 'default_longitude_ssb': parameters.eclipticlongitude, + 'default_latitude_ssb': parameters.eclipticlatitude, + 'default_polarization_ssb': parameters.polarization + } + + def __init__( + self, tc_geo_param=None, longitude_geo_param=None, + latitude_geo_param=None, polarization_geo_param=None, + tc_ssb_param=None, longitude_ssb_param=None, + latitude_ssb_param=None, polarization_ssb_param=None + ): + params = [tc_geo_param, longitude_geo_param, + latitude_geo_param, polarization_geo_param, + tc_ssb_param, longitude_ssb_param, + latitude_ssb_param, polarization_ssb_param] + + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_geo_param = params[0] + self.longitude_geo_param = params[1] + self.latitude_geo_param = params[2] + self.polarization_geo_param = params[3] + self.tc_ssb_param = params[4] + self.longitude_ssb_param = params[5] + self.latitude_ssb_param = params[6] + self.polarization_ssb_param = params[7] + self._inputs = [self.tc_geo_param, self.longitude_geo_param, + self.latitude_geo_param, self.polarization_geo_param] + self._outputs = [self.tc_ssb_param, self.longitude_ssb_param, + self.latitude_ssb_param, self.polarization_ssb_param] + + super(GEOToSSB, self).__init__() + + def transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the geocentric frame to the corresponding + values in the SSB frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_ssb_param], out[self.longitude_ssb_param], \ + out[self.latitude_ssb_param], out[self.polarization_ssb_param] = \ + coordinates.geo_to_ssb( + maps[self.tc_geo_param], maps[self.longitude_geo_param], + maps[self.latitude_geo_param], maps[self.polarization_geo_param] + ) + return self.format_output(maps, out) + + def inverse_transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the SSB frame to the corresponding + values in the geocentric frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_geo_param], out[self.longitude_geo_param], \ + out[self.latitude_geo_param], out[self.polarization_geo_param] = \ + coordinates.ssb_to_geo( + maps[self.tc_ssb_param], maps[self.longitude_ssb_param], + maps[self.latitude_ssb_param], maps[self.polarization_ssb_param] + ) + return self.format_output(maps, out) + + @classmethod + def from_config(cls, cp, section, outputs): + tag = outputs + skip_opts = [] + additional_opts = {} + + # get custom variable names + variables = { + 'tc-geo': cls.default_params_name['default_tc_geo'], + 'longitude-geo': cls.default_params_name['default_longitude_geo'], + 'latitude-geo': cls.default_params_name['default_latitude_geo'], + 'polarization-geo': cls.default_params_name[ + 'default_polarization_geo'], + 'tc-ssb': cls.default_params_name['default_tc_ssb'], + 'longitude-ssb': cls.default_params_name['default_longitude_ssb'], + 'latitude-ssb': cls.default_params_name['default_latitude_ssb'], + 'polarization-ssb': cls.default_params_name[ + 'default_polarization_ssb'] + } + for param_name in variables.keys(): + name_underline = param_name.replace('-', '_') + if cp.has_option("-".join([section, outputs]), param_name): + skip_opts.append(param_name) + additional_opts.update( + {name_underline+'_param': cp.get_opt_tag( + section, param_name, tag)}) + else: + additional_opts.update( + {name_underline+'_param': variables[param_name]}) + + return super(GEOToSSB, cls).from_config( + cp, section, outputs, skip_opts=skip_opts, + additional_opts=additional_opts + ) + + +class LISAToSSB(BaseTransform): + """Converts arrival time, sky localization, and polarization angle in the + LISA frame to the corresponding values in the SSB frame.""" + + name = "lisa_to_ssb" + + default_params_name = { + 'default_tc_lisa': parameters.tc, + 'default_longitude_lisa': parameters.eclipticlongitude, + 'default_latitude_lisa': parameters.eclipticlatitude, + 'default_polarization_lisa': parameters.polarization, + 'default_tc_ssb': parameters.tc, + 'default_longitude_ssb': parameters.eclipticlongitude, + 'default_latitude_ssb': parameters.eclipticlatitude, + 'default_polarization_ssb': parameters.polarization + } + + def __init__( + self, tc_lisa_param=None, longitude_lisa_param=None, + latitude_lisa_param=None, polarization_lisa_param=None, + tc_ssb_param=None, longitude_ssb_param=None, + latitude_ssb_param=None, polarization_ssb_param=None + ): + params = [tc_lisa_param, longitude_lisa_param, + latitude_lisa_param, polarization_lisa_param, + tc_ssb_param, longitude_ssb_param, + latitude_ssb_param, polarization_ssb_param] + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_lisa_param = params[0] + self.longitude_lisa_param = params[1] + self.latitude_lisa_param = params[2] + self.polarization_lisa_param = params[3] + self.tc_ssb_param = params[4] + self.longitude_ssb_param = params[5] + self.latitude_ssb_param = params[6] + self.polarization_ssb_param = params[7] + self._inputs = [self.tc_lisa_param, self.longitude_lisa_param, + self.latitude_lisa_param, self.polarization_lisa_param] + self._outputs = [self.tc_ssb_param, self.longitude_ssb_param, + self.latitude_ssb_param, self.polarization_ssb_param] + super(LISAToSSB, self).__init__() + + def transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the LISA frame to the corresponding + values in the SSB frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_ssb_param], out[self.longitude_ssb_param], \ + out[self.latitude_ssb_param], out[self.polarization_ssb_param] = \ + coordinates.lisa_to_ssb( + maps[self.tc_lisa_param], maps[self.longitude_lisa_param], + maps[self.latitude_lisa_param], maps[self.polarization_lisa_param] + ) + return self.format_output(maps, out) + + def inverse_transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the SSB frame to the corresponding + values in the LISA frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_lisa_param], out[self.longitude_lisa_param], \ + out[self.latitude_lisa_param], \ + out[self.polarization_lisa_param] = \ + coordinates.ssb_to_lisa( + maps[self.tc_ssb_param], maps[self.longitude_ssb_param], + maps[self.latitude_ssb_param], maps[self.polarization_ssb_param] + ) + return self.format_output(maps, out) + + @classmethod + def from_config(cls, cp, section, outputs): + tag = outputs + skip_opts = [] + additional_opts = {} + + # get custom variable names + variables = { + 'tc-lisa': cls.default_params_name['default_tc_lisa'], + 'longitude-lisa': cls.default_params_name[ + 'default_longitude_lisa'], + 'latitude-lisa': cls.default_params_name['default_latitude_lisa'], + 'polarization-lisa': cls.default_params_name[ + 'default_polarization_lisa'], + 'tc-ssb': cls.default_params_name['default_tc_ssb'], + 'longitude-ssb': cls.default_params_name['default_longitude_ssb'], + 'latitude-ssb': cls.default_params_name['default_latitude_ssb'], + 'polarization-ssb': cls.default_params_name[ + 'default_polarization_ssb'] + } + for param_name in variables.keys(): + name_underline = param_name.replace('-', '_') + if cp.has_option("-".join([section, outputs]), param_name): + skip_opts.append(param_name) + additional_opts.update( + {name_underline+'_param': cp.get_opt_tag( + section, param_name, tag)}) + else: + additional_opts.update( + {name_underline+'_param': variables[param_name]}) + + return super(LISAToSSB, cls).from_config( + cp, section, outputs, skip_opts=skip_opts, + additional_opts=additional_opts + ) + + +class LISAToGEO(BaseTransform): + """Converts arrival time, sky localization, and polarization angle in the + LISA frame to the corresponding values in the geocentric frame.""" + + name = "lisa_to_geo" + + default_params_name = { + 'default_tc_lisa': parameters.tc, + 'default_longitude_lisa': parameters.eclipticlongitude, + 'default_latitude_lisa': parameters.eclipticlatitude, + 'default_polarization_lisa': parameters.polarization, + 'default_tc_geo': parameters.tc, + 'default_longitude_geo': parameters.ra, + 'default_latitude_geo': parameters.dec, + 'default_polarization_geo': parameters.polarization + } + + def __init__( + self, tc_lisa_param=None, longitude_lisa_param=None, + latitude_lisa_param=None, polarization_lisa_param=None, + tc_geo_param=None, longitude_geo_param=None, + latitude_geo_param=None, polarization_geo_param=None + ): + params = [tc_lisa_param, longitude_lisa_param, + latitude_lisa_param, polarization_lisa_param, + tc_geo_param, longitude_geo_param, + latitude_geo_param, polarization_geo_param] + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_lisa_param = params[0] + self.longitude_lisa_param = params[1] + self.latitude_lisa_param = params[2] + self.polarization_lisa_param = params[3] + self.tc_geo_param = params[4] + self.longitude_geo_param = params[5] + self.latitude_geo_param = params[6] + self.polarization_geo_param = params[7] + self._inputs = [self.tc_lisa_param, self.longitude_lisa_param, + self.latitude_lisa_param, self.polarization_lisa_param] + self._outputs = [self.tc_geo_param, self.longitude_geo_param, + self.latitude_geo_param, self.polarization_geo_param] + super(LISAToGEO, self).__init__() + + def transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the LISA frame to the corresponding + values in the geocentric frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_geo_param], out[self.longitude_geo_param], \ + out[self.latitude_geo_param], out[self.polarization_geo_param] = \ + coordinates.lisa_to_geo( + maps[self.tc_lisa_param], maps[self.longitude_lisa_param], + maps[self.latitude_lisa_param], maps[self.polarization_lisa_param] + ) + return self.format_output(maps, out) + + def inverse_transform(self, maps): + """This function transforms arrival time, sky localization, + and polarization angle in the geocentric frame to the corresponding + values in the LISA frame. + + Parameters + ---------- + maps : a mapping object + + Returns + ------- + out : dict + A dict with key as parameter name and value as numpy.array or float + of transformed values. + """ + out = {} + out[self.tc_lisa_param], out[self.longitude_lisa_param], \ + out[self.latitude_lisa_param], \ + out[self.polarization_lisa_param] = \ + coordinates.geo_to_lisa( + maps[self.tc_geo_param], maps[self.longitude_geo_param], + maps[self.latitude_geo_param], maps[self.polarization_geo_param] + ) + return self.format_output(maps, out) + + @classmethod + def from_config(cls, cp, section, outputs): + tag = outputs + skip_opts = [] + additional_opts = {} + + # get custom variable names + variables = { + 'tc-lisa': cls.default_params_name['default_tc_lisa'], + 'longitude-lisa': cls.default_params_name[ + 'default_longitude_lisa'], + 'latitude-lisa': cls.default_params_name['default_latitude_lisa'], + 'polarization-lisa': cls.default_params_name[ + 'default_polarization_lisa'], + 'tc-geo': cls.default_params_name['default_tc_geo'], + 'longitude-geo': cls.default_params_name['default_longitude_geo'], + 'latitude-geo': cls.default_params_name['default_latitude_geo'], + 'polarization-geo': cls.default_params_name[ + 'default_polarization_geo'] + } + for param_name in variables.keys(): + name_underline = param_name.replace('-', '_') + if cp.has_option("-".join([section, outputs]), param_name): + skip_opts.append(param_name) + additional_opts.update( + {name_underline+'_param': cp.get_opt_tag( + section, param_name, tag)}) + else: + additional_opts.update( + {name_underline+'_param': variables[param_name]}) + + return super(LISAToGEO, cls).from_config( + cp, section, outputs, skip_opts=skip_opts, + additional_opts=additional_opts + ) + + class Log(BaseTransform): """Applies a log transform from an `inputvar` parameter to an `outputvar` parameter. This is the inverse of the exponent transform. @@ -2060,6 +2456,117 @@ class ChiPToCartesianSpin(CartesianSpinToChiP): inverse_jacobian = inverse.jacobian +class SSBToGEO(GEOToSSB): + """The inverse of GEOToSSB.""" + + name = "ssb_to_geo" + inverse = GEOToSSB + transform = inverse.inverse_transform + inverse_transform = inverse.transform + + def __init__( + self, tc_geo_param=None, longitude_geo_param=None, + latitude_geo_param=None, polarization_geo_param=None, + tc_ssb_param=None, longitude_ssb_param=None, + latitude_ssb_param=None, polarization_ssb_param=None + ): + params = [tc_geo_param, longitude_geo_param, + latitude_geo_param, polarization_geo_param, + tc_ssb_param, longitude_ssb_param, + latitude_ssb_param, polarization_ssb_param] + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_geo_param = params[0] + self.longitude_geo_param = params[1] + self.latitude_geo_param = params[2] + self.polarization_geo_param = params[3] + self.tc_ssb_param = params[4] + self.longitude_ssb_param = params[5] + self.latitude_ssb_param = params[6] + self.polarization_ssb_param = params[7] + self._inputs = [self.tc_ssb_param, self.longitude_ssb_param, + self.latitude_ssb_param, self.polarization_ssb_param] + self._outputs = [self.tc_geo_param, self.longitude_geo_param, + self.latitude_geo_param, self.polarization_geo_param] + + +class SSBToLISA(LISAToSSB): + """The inverse of LISAToSSB.""" + + name = "ssb_to_lisa" + inverse = LISAToSSB + transform = inverse.inverse_transform + inverse_transform = inverse.transform + + def __init__( + self, tc_lisa_param=None, longitude_lisa_param=None, + latitude_lisa_param=None, polarization_lisa_param=None, + tc_ssb_param=None, longitude_ssb_param=None, + latitude_ssb_param=None, polarization_ssb_param=None + ): + params = [tc_lisa_param, longitude_lisa_param, + latitude_lisa_param, polarization_lisa_param, + tc_ssb_param, longitude_ssb_param, + latitude_ssb_param, polarization_ssb_param] + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_lisa_param = params[0] + self.longitude_lisa_param = params[1] + self.latitude_lisa_param = params[2] + self.polarization_lisa_param = params[3] + self.tc_ssb_param = params[4] + self.longitude_ssb_param = params[5] + self.latitude_ssb_param = params[6] + self.polarization_ssb_param = params[7] + self._inputs = [self.tc_ssb_param, self.longitude_ssb_param, + self.latitude_ssb_param, self.polarization_ssb_param] + self._outputs = [self.tc_lisa_param, self.longitude_lisa_param, + self.latitude_lisa_param, self.polarization_lisa_param] + + +class GEOToLISA(LISAToGEO): + """The inverse of LISAToGEO.""" + + name = "geo_to_lisa" + inverse = LISAToGEO + transform = inverse.inverse_transform + inverse_transform = inverse.transform + + def __init__( + self, tc_lisa_param=None, longitude_lisa_param=None, + latitude_lisa_param=None, polarization_lisa_param=None, + tc_geo_param=None, longitude_geo_param=None, + latitude_geo_param=None, polarization_geo_param=None + ): + params = [tc_lisa_param, longitude_lisa_param, + latitude_lisa_param, polarization_lisa_param, + tc_geo_param, longitude_geo_param, + latitude_geo_param, polarization_geo_param] + for index in range(len(params)): + if params[index] is None: + key = list(self.default_params_name.keys())[index] + params[index] = self.default_params_name[key] + + self.tc_lisa_param = params[0] + self.longitude_lisa_param = params[1] + self.latitude_lisa_param = params[2] + self.polarization_lisa_param = params[3] + self.tc_geo_param = params[4] + self.longitude_geo_param = params[5] + self.latitude_geo_param = params[6] + self.polarization_geo_param = params[7] + self._inputs = [self.tc_geo_param, self.longitude_geo_param, + self.latitude_geo_param, self.polarization_geo_param] + self._outputs = [self.tc_lisa_param, self.longitude_lisa_param, + self.latitude_lisa_param, self.polarization_lisa_param] + + class Exponent(Log): """Applies an exponent transform to an `inputvar` parameter. @@ -2202,6 +2709,9 @@ def from_config(cls, cp, section, outputs, ChiPToCartesianSpin.inverse = CartesianSpinToChiP Log.inverse = Exponent Logit.inverse = Logistic +GEOToSSB.inverse = SSBToGEO +LISAToSSB.inverse = SSBToLISA +LISAToGEO.inverse = GEOToLISA # @@ -2241,6 +2751,12 @@ def from_config(cls, cp, section, outputs, LambdaFromTOVFile.name: LambdaFromTOVFile, LambdaFromMultipleTOVFiles.name: LambdaFromMultipleTOVFiles, AlignTotalSpin.name: AlignTotalSpin, + GEOToSSB.name: GEOToSSB, + SSBToGEO.name: SSBToGEO, + LISAToSSB.name: LISAToSSB, + SSBToLISA.name: SSBToLISA, + LISAToGEO.name: LISAToGEO, + GEOToLISA.name: GEOToLISA, } # standard CBC transforms: these are transforms that do not require input @@ -2269,6 +2785,9 @@ def from_config(cls, cp, section, outputs, PrecessionMassSpinToCartesianSpin(), ChiPToCartesianSpin(), ChirpDistanceToDistance(), + GEOToSSB(), + LISAToSSB(), + LISAToGEO(), ] common_cbc_inverse_transforms = [ _t.inverse() diff --git a/test/test_coordinates_space.py b/test/test_coordinates_space.py new file mode 100644 index 00000000000..6e4f8883371 --- /dev/null +++ b/test/test_coordinates_space.py @@ -0,0 +1,358 @@ +# Copyright (C) 2023 Shichao Wu, Alex Nitz +# +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + +import numpy +from pycbc.coordinates import space +import unittest +from utils import simple_exit + + +seed = 8202 +numpy.random.seed(seed) + +# ranges to draw random numbers for each parameter +RANGES = { + "tc_geo" : (1126259462.43, 1426259462.43), + "ra" : (0.0, 2 * numpy.pi), + "dec" : (-numpy.pi / 2, numpy.pi / 2), + "polarization_geo" : (0.0, 2 * numpy.pi), + "tc_ssb" : (1126259462.43, 1426259462.43), + "eclipticlongitude_ssb" : (0.0, 2 * numpy.pi), + "eclipticlatitude_ssb" : (-numpy.pi / 2, numpy.pi / 2), + "polarization_ssb" : (0.0, 2 * numpy.pi), + "tc_lisa" : (1126259462.43, 1426259462.43), + "eclipticlongitude_lisa" : (0.0, 2 * numpy.pi), + "eclipticlatitude_lisa" : (-numpy.pi / 2, numpy.pi / 2), + "polarization_lisa" : (0.0, 2 * numpy.pi), +} + +def almost_equal(derived_val, check_val, precision=1e-2): + """Checks whether the difference in the derived and check values are less + than the given precision. + """ + allpass = numpy.allclose(derived_val, check_val, atol=precision) + if not allpass: + absdiff = abs(derived_val - check_val) + maxidx = absdiff.argmax() + maxdiff = absdiff[maxidx] + else: + maxdiff = maxidx = None + return allpass, maxdiff, maxidx + +def curve_similarity(curve_1, curve_2): + """Using the Euclidean distance to check + the similarity between two curves. + """ + return numpy.linalg.norm(curve_1 - curve_2) + + +class TestParams(unittest.TestCase): + def setUp(self, *args): + self.numtests = 1000 + self.precision = 1e-8 + + # generate some random points + random_params = { + p : numpy.random.uniform(*RANGES[p], size=self.numtests) + for p in RANGES.keys()} + self.tc_ssb = random_params['tc_ssb'] + self.eclipticlongitude_ssb = random_params['eclipticlongitude_ssb'] + self.eclipticlatitude_ssb = random_params['eclipticlatitude_ssb'] + self.polarization_ssb = random_params['polarization_ssb'] + self.tc_lisa = random_params['tc_lisa'] + self.eclipticlongitude_lisa = random_params['eclipticlongitude_lisa'] + self.eclipticlatitude_lisa = random_params['eclipticlatitude_lisa'] + self.polarization_lisa = random_params['polarization_lisa'] + self.tc_geo = random_params['tc_geo'] + self.ra = random_params['ra'] + self.dec = random_params['dec'] + self.polarization_geo = random_params['polarization_geo'] + + # calculate derived parameters from each + + self.tc_lisa_derived, self.eclipticlongitude_lisa_derived, \ + self.eclipticlatitude_lisa_derived, self.polarization_lisa_derived = \ + space.ssb_to_lisa( + self.tc_ssb, self.eclipticlongitude_ssb, + self.eclipticlatitude_ssb, self.polarization_ssb) + + self.tc_geo_derived, self.ra_derived, \ + self.dec_derived, self.polarization_geo_derived = \ + space.lisa_to_geo( + self.tc_lisa, self.eclipticlongitude_lisa, + self.eclipticlatitude_lisa, self.polarization_lisa) + + self.tc_ssb_derived, self.eclipticlongitude_ssb_derived, \ + self.eclipticlatitude_ssb_derived, self.polarization_ssb_derived = \ + space.geo_to_ssb( + self.tc_geo, self.ra, + self.dec, self.polarization_geo) + + def test_round_robin(self): + """Computes inverse transformations to get original parameters from + derived, then compares them to the original. + """ + msg = '{} does not recover same {}; max difference: {}; inputs: {}' + # following lists (function to check, + # arguments to pass to the function, + # name of self's attribute to compare to) + fchecks = [ + (space.ssb_to_geo, + (self.tc_ssb_derived, self.eclipticlongitude_ssb_derived, + self.eclipticlatitude_ssb_derived, + self.polarization_ssb_derived), + ('tc_geo', 'ra', + 'dec', 'polarization_geo')), + (space.geo_to_lisa, + (self.tc_geo_derived, self.ra_derived, + self.dec_derived, self.polarization_geo_derived), + ('tc_lisa', 'eclipticlongitude_lisa', + 'eclipticlatitude_lisa', 'polarization_lisa')), + (space.lisa_to_ssb, + (self.tc_lisa_derived, self.eclipticlongitude_lisa_derived, + self.eclipticlatitude_lisa_derived, + self.polarization_lisa_derived), + ('tc_ssb', 'eclipticlongitude_ssb', + 'eclipticlatitude_ssb', 'polarization_ssb')), + ] + + for func, args, compval in fchecks: + func_tuple = func(*args) + attr_tuple = [getattr(self, i) for i in compval] + for i in range(len(func_tuple)): + derived_val = func_tuple[i] + check_val = attr_tuple[i] + passed, maxdiff, maxidx = almost_equal(derived_val, check_val, + self.precision) + if not passed: + failinputs = [p[maxidx] for p in args] + else: + failinputs = None + self.assertTrue(passed, msg.format( + func, compval, maxdiff, failinputs)) + + def test_polarization(self): + """Set up a custom V-shape ground-based detector which is almost + co-aligned and co-located to LISA-Z channel. When the long-wavelength + approximation is valid, the detector responses of those two should be + very similar. The transform functions in `pycbc.coordinates.space` are + used in the calculation. + """ + from pycbc.waveform import get_fd_det_waveform + from pycbc.coordinates.space import (ssb_to_lisa, lisa_to_ssb, + ssb_to_geo, lisa_to_geo) + + from pycbc.types.frequencyseries import FrequencySeries + from pycbc.detector import (Detector, load_detector_config, + add_detector_on_earth, + _ground_detectors) + from pycbc.waveform import get_td_waveform, get_fd_waveform + import importlib + + def is_module_installed(module_name): + try: + importlib.import_module(module_name) + return True + except ImportError: + return False + + def strain_generator(det='D1', model='IMRPhenomD', fs=4096, df=None, + flow=2, fref=2, tc=0, params=None, + fd_taper=False): + + # Generate a waveform at the detector-frame. + hp, hc = get_td_waveform(approximant=model, + mass1=params['mass1'], mass2=params['mass2'], + spin1x=0, spin1y=0, + spin1z=params['spin1z'], spin2x=0, + spin2y=0, spin2z=params['spin2z'], + distance=params['distance'], + coa_phase=params['coa_phase'], + inclination=params['inclination'], f_lower=flow, + f_ref=fref, delta_t=1.0/fs) + + # Set merger time to 'tc'. + hp.start_time += tc + hc.start_time += tc + + # Project GW waveform onto GW detectors. + ra = params['ra'] + dec = params['dec'] + psi = params['polarization'] + time = hp.start_time + + det_1 = Detector("D1") + fp_1, fc_1 = det_1.antenna_pattern( + right_ascension=ra, declination=dec, + polarization=psi, t_gps=tc) + + # Not take the rotation of the earth into account. + ht_1 = fp_1*hp + fc_1*hc + ht_list = [ht_1] + + return ht_list + + + # Checking if bbhx is installed, if not, then ignore this test + # and raise a warning. + # TODO: we need to implement a better way to install bbhx package. + bbhx_installed = is_module_installed('bbhx') + if not bbhx_installed: + passed = True + print("Ignore polarization test, because bbhx is not installed.") + else: + # cunstom E3 + # All of those hard-coded numbers are carefully chosen, + # in order to almost co-align and co-locate both 2 detectors. + fine_tunning = 7992204.094271309 + OMEGA_0 = 1.99098659277e-7 + yangle = numpy.pi / 2 + fine_tunning * OMEGA_0 + add_detector_on_earth(name='D1', longitude=1.8895427761465164, + latitude=0.11450614784814996, yangle=yangle, + xangle=yangle+numpy.pi/3, height=0) + + # set parameters + params = {} + YRSID_SI = 31558149.763545603 + params['ref_frame'] = 'SSB' + params['approximant'] = 'BBHX_PhenomD' + params['base_approximant'] = 'BBHX_PhenomD' + params['coa_phase'] = 0.0 + params['mass1'] = 1e6 + params['mass2'] = 8e5 + params['spin1z'] = 0.0 + params['spin2z'] = 0.0 + params['distance'] = 410 + params['inclination'] = numpy.pi/2 # edge-on + params['t_obs_start'] = 1*YRSID_SI + params['delta_f'] = 1./params['t_obs_start'] + params['f_lower'] = 1e-4 + params['f_ref'] = 8e-4 + params['f_final'] = 0.1 + params['delta_t'] = 1/0.2 + params['t_offset'] = 9206958.120016199 # 0 degrees + t_lisa = YRSID_SI - params['t_offset'] + fine_tunning + + nx = 100 + longitude_array_high_res = numpy.linspace( + 0, 2*numpy.pi, num=nx, endpoint=False) + + amp_E3_psi = [] + phase_E3_psi = [] + amp_LISA3_psi = [] + phase_LISA3_psi = [] + + for polarization_lisa in numpy.linspace( + 0, 2*numpy.pi, endpoint=False, num=nx): + params['tc'], params['eclipticlongitude'], \ + params['eclipticlatitude'], params['polarization'] = \ + lisa_to_ssb(t_lisa, 0, numpy.pi/4, + polarization_lisa, params['t_offset']) + + nparams = {'mass1':params['mass1'], 'mass2':params['mass2'], + 'spin1z':params['spin1z'], + 'spin2z':params['spin2z'], + 'f_lower':params['f_lower']} + + bbhx_fd = get_fd_det_waveform( + ifos=['LISA_A','LISA_E','LISA_T'], **params) + lisa_X_fd = -numpy.sqrt(2)/2 * bbhx_fd['LISA_A'] + \ + numpy.sqrt(6)/6 * bbhx_fd['LISA_E'] + \ + numpy.sqrt(3)/3 * bbhx_fd['LISA_T'] + lisa_Y_fd = -numpy.sqrt(6)/3 * bbhx_fd['LISA_E'] + \ + numpy.sqrt(3)/3 * bbhx_fd['LISA_T'] + lisa_Z_fd = numpy.sqrt(2)/2 * bbhx_fd['LISA_A'] + \ + numpy.sqrt(6)/6 * bbhx_fd['LISA_E'] + \ + numpy.sqrt(3)/3 * bbhx_fd['LISA_T'] + channel_xyz_fd = {'LISA_X':lisa_X_fd, + 'LISA_Y':lisa_Y_fd, + 'LISA_Z':lisa_Z_fd} + + t_geo, ra, dec, polarization_geo = \ + ssb_to_geo(params['tc'], params['eclipticlongitude'], + params['eclipticlatitude'], + params['polarization']) + + params_3g = params.copy() + params_3g['tc'] = t_geo + params_3g['ra'] = ra + params_3g['dec'] = dec + params_3g['polarization'] = polarization_geo + params_3g['coa_phase'] = numpy.pi/2 - params['coa_phase'] + + E3_signal = strain_generator(det='D1', model='IMRPhenomD', + fs=1./params_3g['delta_t'], + flow=params_3g['f_lower'], + fref=params_3g['f_ref'], + tc=params_3g['tc'], + params=params_3g, + fd_taper=False) + E3_signal_fd = { + 'DIY_E3':E3_signal[0].to_frequencyseries( + params_3g['delta_f']) + } + + index_E3 = numpy.where(numpy.abs( + E3_signal_fd['DIY_E3'].sample_frequencies - + params_3g['f_ref']) < 1e-7) + index_LISA3 = numpy.where(numpy.abs( + bbhx_fd['LISA_A'].sample_frequencies - + params['f_ref']) < 1e-7) + amp_E3 = numpy.abs(E3_signal_fd['DIY_E3'])[index_E3[0][0]] + amp_LISA3 = numpy.abs( + channel_xyz_fd['LISA_Z'][index_LISA3[0][0]]) + phase_E3 = numpy.rad2deg(numpy.angle( + E3_signal_fd['DIY_E3'][index_E3[0][0]])) + phase_LISA3 = numpy.rad2deg(numpy.angle( + channel_xyz_fd['LISA_Z'][index_LISA3[0][0]])) + + amp_E3_psi.append(amp_E3) + phase_E3_psi.append(phase_E3) + amp_LISA3_psi.append(amp_LISA3) + phase_LISA3_psi.append(phase_LISA3) + + dist_amp = curve_similarity(amp_E3_psi/numpy.mean(amp_E3_psi), + amp_LISA3_psi/numpy.mean( + amp_LISA3_psi)) + dist_phase = curve_similarity(numpy.array(phase_E3_psi), + numpy.array(phase_LISA3_psi)) + # Here we compare the amplitude and phase as a function of the + # polarization angle in the LISA frame, because E3 and LISA-Z are + # almost co-aligned and co-located, their detector response should + # be very similar (when long-wavelength approximation holds), + # in our case, `dist_amp` and `dist_phase` should be very similar + # (if the frame transform functions work correctly), users can + # modify this script to plot those curves (amp_E3_psi, + # amp_LISA3_psi, phase_E3_psi, amp_LISA3_psi). The hard-coded + # values below are the reference values for the difference which + # I got on my laptop, those curves are not exactly the same, + # especially for the phase, but they are almost consistent with + # each other visually. + if (numpy.abs(dist_amp - 0.2838670151034317) < 1e-2) and \ + (numpy.abs(dist_phase - 151.77197349820668) < 1e-2): + passed = True + else: + passed = False + + self.assertTrue(passed) + + +suite = unittest.TestSuite() +suite.addTest(unittest.TestLoader().loadTestsFromTestCase(TestParams)) + +if __name__ == '__main__': + results = unittest.TextTestRunner(verbosity=2).run(suite) + simple_exit(results) diff --git a/test/test_transforms.py b/test/test_transforms.py index c5cd0df63e9..27b628cb13a 100644 --- a/test/test_transforms.py +++ b/test/test_transforms.py @@ -44,6 +44,12 @@ "xi1" : (0.0, 1.0), "xi2" : (0.0, 1.0), "chirp_distance" : (2.0, 10.0), + "tc" : (1126259462.43, 1526259462.43), + "ra" : (0.0, 2 * numpy.pi), + "dec" : (-numpy.pi / 2, numpy.pi / 2), + "eclipticlongitude" : (0.0, 2 * numpy.pi), + "eclipticlatitude" : (-numpy.pi / 2, numpy.pi / 2), + "polarization" : (0.0, 2 * numpy.pi), } # tests only need to happen on the CPU @@ -96,4 +102,3 @@ def test_inverse(self): if __name__ == "__main__": results = unittest.TextTestRunner(verbosity=2).run(suite) simple_exit(results) - diff --git a/tox.ini b/tox.ini index e41f67686d2..54317f89982 100644 --- a/tox.ini +++ b/tox.ini @@ -18,12 +18,27 @@ conda_deps= openssl=1.1 m2crypto conda_channels=conda-forge +platform = lin: linux + mac: darwin # This test should run on almost anybody's environment [testenv:py-unittest] deps = {[base]deps} pytest + ; Needed for `BBHx` package to work with PyCBC + git+https://github.com/mikekatz04/BBHx.git; sys_platform == 'linux' + git+https://github.com/gwastro/BBHX-waveform-model.git; sys_platform == 'linux' +conda_deps= + mysqlclient + lin: gcc_linux-64>=12.2.0 + lin: gxx_linux-64>=12.2.0 + ; mac doesn't work, need fix + ; mac: clang_osx-64 + ; mac: clangxx_osx-64 + gsl + lapack==3.6.1 +conda_channels=conda-forge commands = pytest # The following are long running or may require @@ -54,11 +69,14 @@ deps = {[base]deps} ; Needed for `BBHx` package to work with PyCBC git+https://github.com/mikekatz04/BBHx.git; sys_platform == 'linux' - git+https://github.com/ConWea/BBHX-waveform-model.git; sys_platform == 'linux' + git+https://github.com/gwastro/BBHX-waveform-model.git; sys_platform == 'linux' conda_deps= mysqlclient - gcc_linux-64>=12.2.0 - gxx_linux-64>=12.2.0 + lin: gcc_linux-64>=12.2.0 + lin: gxx_linux-64>=12.2.0 + ; mac doesn't work, need fix + ; mac: clang_osx-64 + ; mac: clangxx_osx-64 binutils_linux-64>=2.39 gsl lapack==3.6.1