Coverage for pesummary/gw/conversions/angles.py: 69.4%

1# Licensed under an MIT style license -- see LICENSE.md 

2 

3import numpy as np 

4from pesummary.utils.decorators import array_input 

5 

6__author__ = [ 

7 "Stephen Fairhurst <stephen.fairhurst@ligo.org>", 

8 "Rhys Green <rhys.green@ligo.org>", 

9 "Charlie Hoy <charlie.hoy@ligo.org>", 

10] 

11 

12 

13def _dpsi(theta_jn, phi_jl, beta): 

14 """Calculate the difference between the polarization with respect to the 

15 total angular momentum and the polarization with respect to the orbital 

16 angular momentum 

17 """ 

18 if theta_jn == 0: 

19 return -1. * phi_jl 

20 n = np.array([np.sin(theta_jn), 0, np.cos(theta_jn)]) 

21 j = np.array([0, 0, 1]) 

22 l = np.array([ 

23 np.sin(beta) * np.sin(phi_jl), np.sin(beta) * np.cos(phi_jl), np.cos(beta) 

24 ]) 

25 p_j = np.cross(n, j) 

26 p_j /= np.linalg.norm(p_j) 

27 p_l = np.cross(n, l) 

28 p_l /= np.linalg.norm(p_l) 

29 cosine = np.inner(p_j, p_l) 

30 sine = np.inner(n, np.cross(p_j, p_l)) 

31 dpsi = np.pi / 2 + np.sign(sine) * np.arccos(cosine) 

32 return dpsi 

33 

34 

35@array_input() 

36def _dphi(theta_jn, phi_jl, beta): 

37 """Calculate the difference in the phase angle between J-aligned 

38 and L-aligned frames 

39 

40 Parameters 

41 ---------- 

42 theta_jn: np.ndarray 

43 the angle between J and line of sight 

44 phi_jl: np.ndarray 

45 the precession phase 

46 beta: np.ndarray 

47 the opening angle (angle between J and L) 

48 """ 

49 n = np.column_stack( 

50 [np.repeat([0], len(theta_jn)), np.sin(theta_jn), np.cos(theta_jn)] 

51 ) 

52 l = np.column_stack( 

53 [ 

54 np.sin(beta) * np.cos(phi_jl), np.sin(beta) * np.sin(phi_jl), 

55 np.cos(beta) 

56 ] 

57 ) 

58 cosi = [np.inner(nn, ll) for nn, ll in zip(n, l)] 

59 inc = np.arccos(cosi) 

60 sign = np.sign(np.cos(theta_jn) - (np.cos(beta) * np.cos(inc))) 

61 cos_d = np.cos(phi_jl) * np.sin(theta_jn) / np.sin(inc) 

62 inds = np.logical_or(cos_d < -1, cos_d > 1) 

63 cos_d[inds] = np.sign(cos_d[inds]) * 1. 

64 dphi = -1. * sign * np.arccos(cos_d) 

65 return dphi 

66 

67 

68@array_input() 

69def psi_J(psi_L, theta_jn, phi_jl, beta): 

70 """Calculate the polarization with respect to the total angular momentum 

71 """ 

72 dpsi = [] 

73 for i in range(len(theta_jn)): 

74 dpsi.append(_dpsi(theta_jn[i], phi_jl[i], beta[i])) 

75 psi = psi_L + np.array(dpsi) 

76 return psi