Coverage for pesummary/gw/conversions/spins.py: 63.3%
120 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-09 22:34 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3import numpy as np
4from pesummary.utils.decorators import array_input
6__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
8try:
9 from lalsimulation import (
10 SimInspiralTransformPrecessingWvf2PE,
11 SimInspiralTransformPrecessingNewInitialConditions
12 )
13 from lal import MSUN_SI
14except ImportError:
15 pass
18@array_input()
19def viewing_angle_from_inclination(inclination):
20 """Return the viewing angle of the binary given samples for the source
21 inclination angle. For a precessing system, the source inclination angle
22 is theta_jn: the angle between the total angular momentum J and the line of
23 sight N. For a non-precessing system, the source inclination angle is
24 iota: the angle between the total orbital angular momentum L and the line
25 of sight N
26 """
27 return np.min([inclination, np.pi - inclination], axis=0)
30@array_input()
31def _chi_p(mass1, mass2, S1_perp, S2_perp):
32 """Return chi_p given samples for mass1, mass2, S1_perp, S2_perp
33 """
34 mass_ratio = mass2 / mass1
35 chi_p = np.maximum(
36 S1_perp, (4 * mass_ratio + 3) / (3 * mass_ratio + 4) * mass_ratio
37 * S2_perp
38 )
39 return chi_p
42@array_input()
43def chi_p(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
44 """Return chi_p given samples for mass1, mass2, spin1x, spin1y, spin2x,
45 spin2y
46 """
47 S1_perp = np.linalg.norm([spin1x, spin1y], axis=0)
48 S2_perp = np.linalg.norm([spin2x, spin2y], axis=0)
49 return _chi_p(mass1, mass2, S1_perp, S2_perp)
52@array_input()
53def chi_p_from_tilts(mass1, mass2, a_1, tilt_1, a_2, tilt_2):
54 """Return chi_p given samples for mass1, mass2, a_1, tilt_2, a_2, tilt_2
55 """
56 S1_perp = a_1 * np.sin(tilt_1)
57 S2_perp = a_2 * np.sin(tilt_2)
58 return _chi_p(mass1, mass2, S1_perp, S2_perp)
61@array_input()
62def chi_p_2spin(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
63 """Return the magnitude of a modified chi_p which takes into account
64 precessing spin information from both compact objects given samples for
65 mass1, mass2, spin1x, spin1y, spin2x, spin2y. See Eq.9 of arXiv:2012.02209
66 """
67 chi_p_2spin = np.zeros((2, len(mass1)))
68 S1_perp = mass1**2 * np.array([spin1x, spin1y])
69 S2_perp = mass2**2 * np.array([spin2x, spin2y])
70 S_perp = np.sum([S1_perp, S2_perp], axis=0)
71 S1_perp_mag = np.linalg.norm(S1_perp, axis=0)
72 S2_perp_mag = np.linalg.norm(S2_perp, axis=0)
73 mask = S1_perp_mag >= S2_perp_mag
74 chi_p_2spin[:, mask] = (S_perp / (mass1**2 + S2_perp_mag))[:, mask]
75 chi_p_2spin[:, ~mask] = (S_perp / (mass2**2 + S1_perp_mag))[:, ~mask]
76 return np.linalg.norm(chi_p_2spin, axis=0)
79@array_input()
80def chi_eff(mass1, mass2, spin1z, spin2z):
81 """Return chi_eff given samples for mass1, mass2, spin1z, spin2z
82 """
83 return (spin1z * mass1 + spin2z * mass2) / (mass1 + mass2)
86@array_input()
87def phi_12_from_phi1_phi2(phi1, phi2):
88 """Return the difference in azimuthal angle between S1 and S2 given samples
89 for phi1 and phi2
90 """
91 phi12 = phi2 - phi1
92 if isinstance(phi12, float) and phi12 < 0.:
93 phi12 += 2 * np.pi
94 elif isinstance(phi12, np.ndarray):
95 ind = np.where(phi12 < 0.)
96 phi12[ind] += 2 * np.pi
97 return phi12
100@array_input()
101def phi1_from_spins(spin_1x, spin_1y):
102 """Return phi_1 given samples for spin_1x and spin_1y
103 """
104 phi_1 = np.fmod(2 * np.pi + np.arctan2(spin_1y, spin_1x), 2 * np.pi)
105 return phi_1
108@array_input()
109def phi2_from_spins(spin_2x, spin_2y):
110 """Return phi_2 given samples for spin_2x and spin_2y
111 """
112 return phi1_from_spins(spin_2x, spin_2y)
115@array_input()
116def spin_angles(mass_1, mass_2, inc, spin1x, spin1y, spin1z, spin2x, spin2y,
117 spin2z, f_ref, phase):
118 """Return the spin angles given samples for mass_1, mass_2, inc, spin1x,
119 spin1y, spin1z, spin2x, spin2y, spin2z, f_ref, phase
120 """
121 return_float = False
122 if isinstance(mass_1, (int, float)):
123 return_float = True
124 mass_1 = [mass_1]
125 mass_2 = [mass_2]
126 inc = [inc]
127 spin1x = [spin1x]
128 spin1y = [spin1y]
129 spin1z = [spin1z]
130 spin2x = [spin2x]
131 spin2y = [spin2y]
132 spin2z = [spin2z]
133 f_ref = [f_ref]
134 phase = [phase]
136 data = []
137 for i in range(len(mass_1)):
138 theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2 = \
139 SimInspiralTransformPrecessingWvf2PE(
140 incl=inc[i], m1=mass_1[i], m2=mass_2[i], S1x=spin1x[i],
141 S1y=spin1y[i], S1z=spin1z[i], S2x=spin2x[i], S2y=spin2y[i],
142 S2z=spin2z[i], fRef=float(f_ref[i]), phiRef=float(phase[i]))
143 data.append([theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2])
144 if return_float:
145 return data[0]
146 return data
149@array_input()
150def component_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
151 mass_2, f_ref, phase):
152 """Return the component spins given samples for theta_jn, phi_jl, tilt_1,
153 tilt_2, phi_12, a_1, a_2, mass_1, mass_2, f_ref, phase
154 """
155 data = []
156 for i in range(len(theta_jn)):
157 iota, S1x, S1y, S1z, S2x, S2y, S2z = \
158 SimInspiralTransformPrecessingNewInitialConditions(
159 theta_jn[i], phi_jl[i], tilt_1[i], tilt_2[i], phi_12[i],
160 a_1[i], a_2[i], mass_1[i] * MSUN_SI, mass_2[i] * MSUN_SI,
161 float(f_ref[i]), float(phase[i]))
162 data.append([iota, S1x, S1y, S1z, S2x, S2y, S2z])
163 return data
166@array_input()
167def spin_angles_from_azimuthal_and_polar_angles(
168 a_1, a_2, a_1_azimuthal, a_1_polar, a_2_azimuthal, a_2_polar):
169 """Return the spin angles given samples for a_1, a_2, a_1_azimuthal,
170 a_1_polar, a_2_azimuthal, a_2_polar
171 """
172 spin1x = a_1 * np.sin(a_1_polar) * np.cos(a_1_azimuthal)
173 spin1y = a_1 * np.sin(a_1_polar) * np.sin(a_1_azimuthal)
174 spin1z = a_1 * np.cos(a_1_polar)
176 spin2x = a_2 * np.sin(a_2_polar) * np.cos(a_2_azimuthal)
177 spin2y = a_2 * np.sin(a_2_polar) * np.sin(a_2_azimuthal)
178 spin2z = a_2 * np.cos(a_2_polar)
180 data = [[s1x, s1y, s1z, s2x, s2y, s2z] for s1x, s1y, s1z, s2x, s2y, s2z in
181 zip(spin1x, spin1y, spin1z, spin2x, spin2y, spin2z)]
182 return data
185@array_input()
186def tilt_angles_and_phi_12_from_spin_vectors_and_L(a_1, a_2, Ln):
187 """Return the tilt angles and phi_12 given samples for the spin vectors
188 and the orbital angular momentum
190 Parameters
191 ----------
192 a_1: np.ndarray
193 Spin vector for the larger object
194 a_2: np.ndarray
195 Spin vector for the smaller object
196 Ln: np.ndarray
197 Orbital angular momentum of the binary
198 """
199 a_1_norm = np.linalg.norm(a_1)
200 a_2_norm = np.linalg.norm(a_2)
201 Ln /= np.linalg.norm(Ln)
202 a_1_dot = np.dot(a_1, Ln)
203 a_2_dot = np.dot(a_2, Ln)
204 a_1_perp = a_1 - a_1_dot * Ln
205 a_2_perp = a_2 - a_2_dot * Ln
206 cos_tilt_1 = a_1_dot / a_1_norm
207 cos_tilt_2 = a_2_dot / a_2_norm
208 cos_phi_12 = np.dot(a_1_perp, a_2_perp) / (
209 np.linalg.norm(a_1_perp) * np.linalg.norm(a_2_perp)
210 )
211 # set quadrant of phi12
212 phi_12 = np.arccos(cos_phi_12)
213 if np.sign(np.dot(Ln, np.cross(a_1, a_2))) < 0.:
214 phi_12 = 2. * np.pi - phi_12
216 return np.arccos(cos_tilt_1), np.arccos(cos_tilt_2), phi_12
219@array_input()
220def opening_angle(
221 mass_1, mass_2, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, f_ref, phase
222):
223 """Return the opening angle of the system given samples for mass_1, mass_2,
224 cartesian spins and a reference frequency
225 """
226 data = []
227 for i in range(len(mass_1)):
228 beta, _, _, _, _, _, _ = \
229 SimInspiralTransformPrecessingNewInitialConditions(
230 0., phi_jl[i], tilt_1[i], tilt_2[i], phi_12[i],
231 a_1[i], a_2[i], mass_1[i] * MSUN_SI, mass_2[i] * MSUN_SI,
232 float(f_ref[i]), float(phase[i]))
233 data.append(beta)
234 return data