Commit 06711516 authored by Sylvia Biscoveanu's avatar Sylvia Biscoveanu

Adding time delay

Merge branch 'master' of git.ligo.org:Monash/peyote
parents 8aacd8a3 b55d1cab
import numpy as np
import pdb
class likelihood:
def __init__(self, Interferometers, source) :
self.Interferometers = Interferometers
self.source = source
def logL(self, params) :
logL = 0
waveform_polarizations = self.source.model(params)
for Interferometer in self.Interferometers :
class likelihood:
def __init__(self, Interferometers, source):
self.Interferometers = Interferometers
self.source = source
for mode in params['modes'] :
def logL(self, times, params):
logL = 0
waveform_polarizations = self.source.model(times, params)
for Interferometer in self.Interferometers:
for mode in params['modes']:
#det_response = Interferometer.antenna_response(
# params['ra'], params['dec'],
# params['geocent_time'], params['psi'],
# mode )
det_response = 1
det_response = Interferometer.antenna_response(
params['ra'], params['dec'],
params['geocent_time'], params['psi'],
mode)
waveform_polarizations[mode] *= det_response
waveform_polarizations[mode] *= det_response
signal_IFO = np.sum( waveform_polarizations.values() )
signal_IFO = np.sum(waveform_polarizations.values())
#time_shift = Interferometer.time_shift(source.params['geocent_time'])
#signal *= np.exp(-1j*2*np.pi*time_shift) # This is just here as a reminder that a tc shift needs to be performed
# on frequency-domain GWs
#time_shift = Interferometer.time_shift(source.params['geocent_time'])
#signal *= np.exp(-1j*2*np.pi*time_shift)
# This is just here as a reminder that a tc shift needs to be performed
# on frequency-domain GWs
logL += 4. * params['deltaF'] * np.vdot( Interferometer.data - signal_IFO, ( Interferometer.data - signal_IFO ) / Interferometer.psd )
logL += 4. * params['deltaF'] * np.vdot(
Interferometer.data - signal_IFO,
(Interferometer.data - signal_IFO) / Interferometer.psd)
return logL
return logL
import peyote
import numpy as np
import os.path
from astropy.table import Table
from peyote.utils import sampling_frequency, nfft
try:
import lal
import lalsimulation as lalsim
except ImportError:
print("lal is not installed")
class Source:
def __init__(self, name):
......@@ -19,9 +26,48 @@ class SimpleSinusoidSource(Source):
"""
def model(self, parameters):
return {'+': parameters['A'] * np.sin(
parameters['f'] * parameters['geocent_time'])}
def model(self, times, parameters):
return {'plus': parameters['A'] * (
np.sin(parameters['f'] * times)
+ 1j * np.cos(parameters['f'] * times))}
class BinaryBlackHole(Source):
"""
A way of getting a BBH waveform from lal
"""
def waveform_from_lal(self, frequencies, parameters):
luminosity_distance = parameters['luminosity_distance'] * 1e6 * lal.PC_SI
mass_1 = parameters['mass_1'] * lal.MSUN_SI
mass_2 = parameters['mass_2'] * lal.MSUN_SI
longitude_ascending_nodes = 0.0
eccentricity = 0.0
meanPerAno = 0.0
waveform_dictionary = lal.CreateDict()
approximant = lalsim.GetApproximantFromString(parameters['waveform_approximant'])
# delta_f = frequencies[1]-frequencies[0]
# minimum_frequency = np.min(frequencies)
# print(minimum_frequency)
# maximum_frequency = np.max(frequencies)
hplus,hcross = lalsim.SimInspiralChooseFDWaveform(mass_1, mass_2,
parameters['spin_1'][0], parameters['spin_1'][1], parameters['spin_1'][2],
parameters['spin_2'][0], parameters['spin_2'][1], parameters['spin_2'][2],
luminosity_distance, parameters['inclination_angle'],
parameters['waveform_phase'],
longitude_ascending_nodes, eccentricity, meanPerAno,
parameters['deltaF'], parameters['fmin'], parameters['fmax'],
parameters['reference_frequency'],
waveform_dictionary, approximant)
h_plus = hplus.data.data
h_cross = hcross.data.data
return h_plus, h_cross
class Glitch(Source):
......@@ -58,12 +104,13 @@ class Supernova(AstrophysicalSource):
AstrophysicalSource.__init__(self, name, right_ascension, declination, luminosity_distance)
class BinaryBlackHole(CompactBinaryCoalescence):
def __init__(self, name, right_ascension, declination, luminosity_distance, mass_1, mass_2, spin_1, spin_2,
coalescence_time, inclination_angle, waveform_phase, polarisation_angle, eccentricity):
CompactBinaryCoalescence.__init__(self, name, right_ascension, declination, luminosity_distance, mass_1,
mass_2, spin_1, spin_2, coalescence_time, inclination_angle, waveform_phase,
polarisation_angle, eccentricity)
# class BinaryBlackHole(CompactBinaryCoalescence):
# def __init__(self, name, right_ascension, declination, luminosity_distance, mass_1, mass_2, spin_1, spin_2,
# coalescence_time, inclination_angle, waveform_phase, polarisation_angle, eccentricity):
# CompactBinaryCoalescence.__init__(self, name, right_ascension, declination, luminosity_distance, mass_1,
# mass_2, spin_1, spin_2, coalescence_time, inclination_angle, waveform_phase,
# polarisation_angle, eccentricity)
class BinaryNeutronStar(CompactBinaryCoalescence):
......@@ -92,19 +139,24 @@ class BinaryNeutronStarMergerNumericalRelativity(Source):
takes parameters mean_mass, mass_ratio and equation_of_state, directory_path
model takes one parameter `parameters`, a dictionary of Parameters and
returns the waveform model.
returns time,hplus,hcross,freq,Hplus(freq),Hcross(freq)
"""
def model(self, parameters):
mean_mass_string = '{:.0f}'.format(parameters['mean_mass'].value * 1000)
eos_string = parameters['equation_of_state'].value
mass_ratio_string = '{:.0f}'.format(parameters['mass_ratio'].value * 10)
directory_path = parameters['directory_path'].value
mean_mass_string = '{:.0f}'.format(parameters['mean_mass'] * 1000)
eos_string = parameters['equation_of_state']
mass_ratio_string = '{:.0f}'.format(parameters['mass_ratio'] * 10)
directory_path = parameters['directory_path']
file_name = '{}-q{}-M{}.csv'.format(eos_string, mass_ratio_string, mean_mass_string)
full_filename = '{}/{}'.format(directory_path, file_name)
print(full_filename)
if os.path.isfile(file_name):
return full_filename
if not os.path.isfile(full_filename):
print('{} does not exist'.format(full_filename)) # add exception
return(-1)
else: # ok file exists
strain_table = Table.read(full_filename)
Hplus, _ = nfft(strain_table["hplus"], sampling_frequency(strain_table['time']))
Hcross, frequency = nfft(strain_table["hcross"], sampling_frequency(strain_table['time']))
return(strain_table['time'],strain_table["hplus"],strain_table["hcross"],frequency,Hplus,Hcross)
This diff is collapsed.
......@@ -18,17 +18,18 @@ time = utils.create_time_series(sampling_frequency, time_duration)
signal_amplitude = 1e-21
signal_frequency = 100
params = dict(A=signal_amplitude, f=2.*np.pi*signal_frequency, geocent_time=time,
modes=['+'], ra=1, dec=2, psi=0, deltaF=sampling_frequency)
params = dict(A=signal_amplitude, f=2.*np.pi*signal_frequency, geocent_time=1,
modes=['plus'], ra=1, dec=2, psi=0, deltaF=sampling_frequency)
foo = src.SimpleSinusoidSource('foo')
ht_signal = foo.model(params)['+']
ht_signal = foo.model(time, params)['plus']
hf_signal, ff = utils.nfft(ht_signal, sampling_frequency)
"""
Create a noise realisation with a default power spectral density
"""
PSD = det.PowerSpectralDensity() # instantiate a detector psd
PSD.import_power_spectral_density() # import default psd
#PSD.import_power_spectral_density(spectral_density_file="CE_psd.txt") # import cosmic explorer
......@@ -48,8 +49,10 @@ plt.legend(loc='best')
plt.tight_layout()
plt.show()
hf_noise[0] = np.max(hf_noise)
hf_noise[-1] = np.max(hf_noise)
IFO = peyote.detector.H1
IFO.data = hf_signal
IFO.psd = 1
IFO.psd = hf_noise
likelihood = peyote.likelihood.likelihood([IFO], foo)
print likelihood.logL(params)
print(likelihood.logL(time, params))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment