Added functions to compute strain from modes

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}'

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)

watpy/wave/wave.py

 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