diff --git a/minke/sources.py b/minke/sources.py index 1c9fac5623410f8b7cb6a537822523a6fd4a2c92..e45cc52c05ef211ee50bacd92535594d22eca777 100644 --- a/minke/sources.py +++ b/minke/sources.py @@ -219,11 +219,11 @@ class Waveform(object): return row - def interpolate(self, x_old, y_old, x_new): + def interpolate(self, x_old, y_old, x_new, method="linear"): """ Convenience funtion to avoid repeated code """ - interpolator = interp.interp1d(x_old, y_old) + interpolator = interp.interp1d(x_old, y_old, method) return interpolator(x_new) @@ -1073,96 +1073,6 @@ class Yakunin10(Supernova): -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): """ @@ -1217,11 +1127,23 @@ class ADI(LongDuration): 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) + decomposed = self.decompose(filepath) np.savetxt(decomposed_path, decomposed, header="time\thplus\thcross", fmt='%.8e') + #decomposed_path = filepath self.params['phi']=0 self.params['incl']=90 self.params['numrel_data'] = decomposed_path + + def _generate(self, rate = 16384.0, half=False): + data = np.genfromtxt(self.params['numrel_data']) + nsamp = len(data) + hp = lal.CreateREAL8TimeSeries("inj time series", lal.LIGOTimeGPS(0,0), 0, 1.0/rate, lal.StrainUnit, nsamp) + hx = lal.CreateREAL8TimeSeries("inj time series", lal.LIGOTimeGPS(0,0), 0, 1.0/rate, lal.StrainUnit, nsamp) + + hp.data.data = data[:,1] + hx.data.data = data[:,2] + + return hp, hx, np.copy(hp), np.copy(hx) def decompose(self, numrel_file, sample_rate = 16384.0, step_back = 0.01): """ @@ -1268,13 +1190,15 @@ class ADI(LongDuration): end = len(data['hp']) * 1.0 / fs # Make the time array times = np.arange(start, end, 1.0/fs) - - # Load the hp components - strainp = data['hp'] - strainc = data['hc'] # Make the new time vector for the requried sample rate target_times = np.arange(times[0], times[-1], 1.0/sample_rate) + #print len(target_times) + # Load the hp components + strainp = data['hp'].T[0].astype(np.float32) + strainc = data['hc'].T[0].astype(np.float32) + #del data + # Prepare the output matrix output = np.zeros((len(target_times), 3)) diff --git a/tests/data/adi_test.mat b/tests/data/adi_test.mat new file mode 100644 index 0000000000000000000000000000000000000000..0f48b2a39090e1b32e58b190661d97b993211a75 Binary files /dev/null and b/tests/data/adi_test.mat differ diff --git a/tests/test_sources.py b/tests/test_sources.py index fad42405be144193b804bc4c0544afa756d03c27..298aee66d142d95645b3c3fc063d5817128ac5e6 100755 --- a/tests/test_sources.py +++ b/tests/test_sources.py @@ -129,6 +129,40 @@ class TestMinkeAdHocSources(unittest.TestCase): np.testing.assert_array_almost_equal(data[0].data.data[::5000], sgdata) +class TestMinkeADISources(unittest.TestCase): + def setUp(self): + """ + Set things up for the tests by making the MDC set, and defining the various parameter distributions. + """ + + self.datafiles = {} + self.datafiles['ADI'] = 'tests/data/adi_test.mat' + + mdctools.mkdir("./testout/frames") + + self.mdcset = mdctools.MDCSet(['L1', 'H1']) + self.times = distribution.even_time(start = 1126620016, stop = 1136995216, rate = 630720, jitter = 20) + self.angles = distribution.supernova_angle(len(self.times)) + + def test_ADIWaveform(self): + """Check the generation of the ADI waveform.""" + + adi = sources.ADI(time = 100, filepath=self.datafiles['ADI']) + adiwave = adi._generate()[0].data.data[::20000] + reg = np.array([ 0.00000000e+00, 6.62293373e-22, 1.52831641e-21, + 1.54044027e-21, 1.27707122e-21, -2.19617692e-22, + -1.52414321e-21, 4.69396022e-22, 1.21479442e-21, + -1.26396697e-21, -5.87290821e-22, 1.36010407e-21, + 1.18247173e-21, 3.76119821e-22, 3.08731444e-22, + 1.14768993e-21, 1.08488237e-21, -1.50088643e-21, + 1.45463205e-21, 1.67209738e-23, -8.44421117e-22, + 9.54960410e-23, 1.29069231e-21, -9.77430178e-22, + -1.34761100e-21, -1.47183099e-21, 4.61589317e-22, + -1.18249701e-21, -3.27337305e-22, -1.45357383e-21, + 1.45374218e-21, 4.60667061e-22]) + + np.testing.assert_array_almost_equal(adiwave, reg) + class TestMinkeSupernovaSources(unittest.TestCase): def setUp(self): """