From 1a2dc2b50b98c018361834f34b3c0f47b43980eb Mon Sep 17 00:00:00 2001
From: Daniel Williams <daniel.williams@ligo.org>
Date: Thu, 8 Mar 2018 15:08:26 +0000
Subject: [PATCH] Added a new Yakunin class to address #10

---
 minke/sources.py | 92 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 91 insertions(+), 1 deletion(-)

diff --git a/minke/sources.py b/minke/sources.py
index 5bb6171..9e05ef0 100644
--- a/minke/sources.py
+++ b/minke/sources.py
@@ -992,7 +992,6 @@ class Dimmelmeier08(Supernova):
         return output
 
 
-
 class Ringdown(Waveform):
     """
     A class to handle Rindown waveforms.
@@ -1001,7 +1000,98 @@ class Ringdown(Waveform):
     waveform = "GenericRingdown"
 
         
+class Yakunin10(Supernova):
+    """
+    The Yakunin10 waveform.
+    """
+
+    waveform = "Yakunin+10"
+
+    def __init__(self, time, sky_dist=uniform_sky, filepath="s12-time-rhplus_matter.txt", decomposed_path=None, ):
+        """
+
+        Parameters
+        ----------
+        time : float or list 
+           The time period over which the injection should be made. If
+           a list is given they should be the start and end times, and
+           the waveform will be produced at some random point in that
+           time range. If a float is given then the injection will be
+           made at that specific time.
+
+        sky_dist : func
+           The function describing the sky distribution which the injections
+           should be made over. Defaults to a uniform sky.
+
+        filepath : str
+           The filepath to the numerical relativity waveform.
+
+        decomposed_path : str
+           The location where the decomposed waveform file should be stored. Optional.
+        """
+        
+        self._clear_params()
+        self.time = time
+        self.sky_dist = sky_dist
+        if not decomposed_path : decomposed_path = filepath+".dec"
+        if not os.path.isfile(decomposed_path) :
+            decomposed = self.decompose(filepath, sample_rate = 16384.0, step_back = 0.01)
+            np.savetxt(decomposed_path, decomposed, header="time (2,-2) (2,-1) (2,0) (2,1) (2,2)", fmt='%.8e')
+        self.params['phi']=0
+        self.params['incl']=90
+        self.params['numrel_data'] = decomposed_path
+        
+    def decompose(self, numrel_file, sample_rate = 16384.0, step_back = 0.01):
+        """
+        Produce the spherial harmonic decompositions of the Dimmelmeier numerical
+        waveform. This is a special case since it is axisymmetric.
+        
+        Parameters
+        ----------
+        numrel_file : str
+           The location of the numerical relativity waveform file.
+        
+        sample_rate : float
+           The sample rate of the NR file. Defaults to 16384.0 Hz.
+        
+        step_back : float
+           The amount of time, in seconds, of the data which should be included
+           before the peak amplitude. Defaults to 0.01 sec.
+
+        Returns
+        -------
+        decomposition : ndarray
+           The l=2 mode spherical decompositions of the waveform. 
+        """
+        extract_dist = 10e-3
+        # Load the times from the file
+        data = np.loadtxt(numrel_file)
+        data = data.T
+        times = data[0]
+        times -= times[0]
+        
+        # Load the hp components   
+        strain = data[1]
+        # Make the new time vector for the requried sample rate
+        target_times = np.arange(times[0], times[-1], 1.0/sample_rate)
+
+        # Prepare the output matrix
+        output = np.zeros((len(target_times), 11))
+
+        # Add the times in to the first column of said matrix
+        output[:, 0] = target_times #/ lal.MTSUN_SI
+        #
+        # Resample to uniform spacing at 16384 kHz
+        #
+        strain_new = self.interpolate(times, strain, target_times)
+        #
+        # Make the output, and rescale it into dimensionless strain values
+        #
+        output[:,5] = strain_new #/*  ( extract_dist * lal.PC_SI * 1.0e6) 
 
+        return output
+
+    
 class LongDuration(Supernova):
     """
 
-- 
GitLab