From bbdacd0c3d5a5e86817ac73582c1610e103dc985 Mon Sep 17 00:00:00 2001
From: alejandraglz <alejandra.gonzalez@uni-jena.de>
Date: Tue, 1 Nov 2022 12:24:24 +0100
Subject: [PATCH] Added functions to compute strain from modes
---
watpy/utils/units.py | 4 ++++
watpy/wave/gwutils.py | 41 ++++++++++++++++++++++++++++++++++-------
watpy/wave/wave.py | 37 ++++++++++++++++++++++++++++++++++---
3 files changed, 72 insertions(+), 10 deletions(-)
diff --git a/watpy/utils/units.py b/watpy/utils/units.py
index 48d2244..d6af39c 100644
--- a/watpy/utils/units.py
+++ b/watpy/utils/units.py
@@ -5,11 +5,15 @@
uts= {
'Msun_sec': 4.925794970773135e-06,
+ 'Msun_meter': 1.476625061404649406193430731479084713e3
}
def MSun_sec():
return uts['Msun_sec']
+def MSun_meter():
+ return uts['Msun_meter']
+
# For plots:
TEX_KM = r'\mathrm{km}'
diff --git a/watpy/wave/gwutils.py b/watpy/wave/gwutils.py
index bcd3f4b..42e5ee1 100644
--- a/watpy/wave/gwutils.py
+++ b/watpy/wave/gwutils.py
@@ -3,6 +3,7 @@ import scipy as sp
import math
from ..utils import num as num
import warnings as wrn
+from scipy.special import factorial as fact
try:
from scipy import factorial2
except:
@@ -13,6 +14,32 @@ except:
# Waveform analysis utilities
# ------------------------------------------------------------------
+def wigner_d_function(l, m, s, incl):
+ """
+ Wigner-d functions
+ """
+ costheta = np.cos(incl*0.5)
+ sintheta = np.sin(incl*0.5)
+ norm = np.sqrt( (fact(l+m) * fact(l-m) * fact(l+s) * fact(l-s)) )
+ ki = np.amax([0,m-s])
+ kf = np.amin([l+m,l-s])
+ k = np.arange(ki,kf+1)
+ div = 1.0/( fact(k) * fact(l+m-k) * fact(l-s-k) * fact(s-m+k) );
+ dWig = div*( np.power(-1.,k) * np.power(costheta,2*l+m-s-2*k) * np.power(sintheta,2*k+s-m) )
+ return norm * np.sum(dWig)
+
+def spinw_spherical_harm(s, l, m, incl,phi):
+ """
+ Spin-weighted spherical harmonics
+ E.g. https://arxiv.org/abs/0709.0093
+ """
+ if ((l<0) or (m<-l) or (m>l)):
+ raise ValueError("wrong (l,m)")
+ c = np.power(-1.,-s) * np.sqrt( (2.*l+1.)/(4.*np.pi) )
+ dWigner = c * wigner_d_function(l,m,-s,incl)
+ exp_mphi = ( np.cos(m*phi) + 1j * np.sin(m*phi) )
+ return dWigner * exp_mphi
+
def interp_fd_wave(fnew, f, h, kind = 'quadratic'):
"""
Interpolate frequency-domain waveform
@@ -295,7 +322,7 @@ def richardson_extrap_series(p, y, t, h):
ye = np.zeros(n); err = np.zeros(n)
# interpolate all data sets on time grid of extrapolated result
- intrp_y = [ num.linterp(te, t[k], y[k]) for k in range(N-1) ]
+ intrp_y = [ np.interp(te, t[k], y[k]) for k in range(N-1) ]
intrp_y.append(y[-1])
# compute Richardson extrapolation pointwise
@@ -445,20 +472,20 @@ def waveform2energetics(h, h_dot, t, modes, mmodes):
* modes : (l,m) indexes
* mmodes : m indexes
"""
- print('got in')
+
dtime = t[1]-t[0]
- print('step 1')
+
E_GW_dot_ff = {}
E_GW_ff = {}
J_GW_dot_ff = {}
J_GW_ff = {}
for l, m in modes:
- print('step 2')
+
E_GW_dot_ff[(l,m)] = 1.0/(16.*np.pi) * np.abs(h_dot[l,m])**2
E_GW_ff[(l,m)] = num.integrate(E_GW_dot_ff[l,m]) * dtime
J_GW_dot_ff[(l,m)] = 1.0/(16.*np.pi) * m * np.imag(h[l,m] * np.conj(h_dot[l,m]))
J_GW_ff[(l,m)] = num.integrate(J_GW_dot_ff[l,m]) * dtime
- print('step 3')
+
E_GW_dot = {}
E_GW = {}
J_GW_dot = {}
@@ -468,13 +495,13 @@ def waveform2energetics(h, h_dot, t, modes, mmodes):
E_GW[m] = np.zeros_like(t)
J_GW_dot[m] = np.zeros_like(t)
J_GW[m] = np.zeros_like(t)
- print('step 4')
+
for l, m in modes:
E_GW_dot[m] += mnfactor(m) * E_GW_dot_ff[l,m]
E_GW[m] += mnfactor(m) * E_GW_ff[l,m]
J_GW_dot[m] += mnfactor(m) * J_GW_dot_ff[l,m]
J_GW[m] += mnfactor(m) * J_GW_ff[l,m]
- print('step 5')
+
E_GW_dot_all = np.zeros_like(t)
E_GW_all = np.zeros_like(t)
J_GW_dot_all = np.zeros_like(t)
diff --git a/watpy/wave/wave.py b/watpy/wave/wave.py
index 9ad911d..041d8c9 100644
--- a/watpy/wave/wave.py
+++ b/watpy/wave/wave.py
@@ -1,7 +1,8 @@
from ..utils.ioutils import *
from ..utils.num import diff1, diffo
from ..utils.viz import wplot
-from .gwutils import fixed_freq_int_2, waveform2energetics, ret_time
+from ..utils.units import *
+from .gwutils import fixed_freq_int_2, waveform2energetics, ret_time, spinw_spherical_harm
import numpy as np
@@ -353,6 +354,37 @@ class mwaves(object):
fname = "EJ_r"+rad_str+".txt"
np.savetxt('{}/{}'.format(path_out,fname), data, header=headstr)
+ def hlm_to_strain(self,phi=0,inclination=0,add_negative_modes=False):
+ """
+ Build strain from time-domain modes in mass rescaled, geom. units
+ Return result in SI units
+
+ Negative m-modes can be added using positive m-modes,
+ h_{l-m} = (-)^l h^{*}_{lm}
+ """
+ PC_SI = 3.085677581491367e+16 # m
+ MPC_SI = 1e6 * PC_SI
+ wave22 = self.get(l=2, m=2)
+ time = wave22.time_ret() * MSun_sec()
+ distance = wave22.prop['detector.radius']*MPC_SI
+ amplitude_prefactor = wave22.prop['mass'] * MSun_meter() / distance
+ h = np.zeros_like( 1j* time )
+ for (l,m) in self.modes:
+ wavelm = self.get(l=l,m=m)
+ amplitude = wavelm.amplitude(var='h')
+ philm = wavelm.phase(var='h')
+ sYlm = spinw_spherical_harm(-2, l, m, phi, inclination)
+ Alm = amplitude_prefactor * amplitude
+ hlm = Alm * np.exp( - 1j * philm )
+ h += hlm * sYlm
+ if (add_negative_modes):
+ sYlm_neg = spinw_spherical_harm(-2, l, -m, phi, inclination)
+ hlm_neg = (-1)**l * np.conj(hlm)
+ h += hlm_neg * sYlm_neg
+ hplus = np.real(h)
+ hcross = - np.imag(h)
+ return time, hplus, hcross
+
class wave(object):
@@ -622,5 +654,4 @@ class wave(object):
dt = self.time[1] - self.time[0]
return win * fixed_freq_int_2( win * self.p4, fcut, dt = dt)
else:
- return self.h
-
+ return self.h
\ No newline at end of file
--
GitLab