Commit ea68d282 authored by Marc Arene's avatar Marc Arene
Browse files

sigpar_utils.py removed

parent 9fb7c32d
import numpy as np
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
import Library.BNS_Utility_Functions as buf
import time
import bilby.core.utils as bc_utils
import bilby.gw.utils as gw_utils
import Library.config as conf
import Library.set_injection_parameters as set_inj
def binaryFile_to_sigpar(binary_file_name):
"""
Reads a binary.in file and returns the sigpar array and n
"""
params = np.loadtxt(binary_file_name)
(RA, Dec, inc, psi, phic, dist, m1, m2, GMST, f_low, Rmin, m1min, m1max, m2min, m2max, Dmin, Dmax)=params
sigpar, n = buf.CBC_Parameter_Allocation(m1, m2, Dec, RA, inc, psi, phic, dist, GMST, Rmin, m1min+m2min, f_low)
return sigpar, n
def sigpar_to_injection_parameters(sigpar, start_time):
"""
From the hmc sigpar array of parameters, returns the injection_parameters needed by bilby to run.
"""
# bilby.gw.utils.time_delay_geocentric(H.vertex, np.array([0,0,0]), injection_parameters['ra'], injection_parameters['dec'], start_time)
# time_shift_H1 = 0.012139641777406887
time_shift_H1 = gw_utils.time_delay_geocentric(conf.H1_vertex, np.array([0,0,0]), sigpar[7], sigpar[6], start_time)
F_plus_H1 = conf.H1.antenna_response(sigpar[7], sigpar[6], start_time + sigpar[8] - time_shift_H1, sigpar[2], 'plus')
F_cross_H1 = conf.H1.antenna_response(sigpar[7], sigpar[6], start_time + sigpar[8] - time_shift_H1, sigpar[2], 'cross')
phase_shift_H1 = np.angle((0.5*(1+np.cos(sigpar[0])**2)*F_plus_H1) - 1j * np.cos(sigpar[0])*F_cross_H1)
# if phase_shift_H1 < 0: phase_shift_H1 += 2 * np.pi
injection_parameters = dict(
mass_1=sigpar[9],
mass_2=sigpar[10],
chi_1=0.0,
chi_2=0.0,
luminosity_distance=sigpar[3]/MPC,
iota=sigpar[0],
theta_jn=sigpar[0],
psi=sigpar[2],
# phase=0,
# phase=0.5 * (sigpar[1] - 2*np.pi*time_shift_H1/(sigpar[8] - time_shift_H1)),
phase=0.5 * (sigpar[1]),
# geocent_time=geocent_time,
geocent_time=start_time + sigpar[8] - time_shift_H1,#geocent_time,
# geocent_time=start_time + sigpar[8],#geocent_time,
ra=sigpar[7],
dec=sigpar[6],
lambda_1=0.,
lambda_2=0.
)
# print("\n2 * phase = {} * pi".format(2 * injection_parameters['phase']/np.pi))
# import IPython; IPython.embed();
return injection_parameters
def injection_parameters_to_sigpar(injection_parameters, minimum_frequency):
# time_shift_H1 = gw_utils.time_delay_geocentric(conf.H1_vertex, np.array([0,0,0]), sigpar[7], sigpar[6], start_time)
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 180 / np.pi
if 'iota' in injection_parameters.keys():
inc = injection_parameters['iota'] * 180 / np.pi
else:
inc = injection_parameters['theta_jn'] * 180 / np.pi
psi = injection_parameters['psi'] * 180 / np.pi
phic = (injection_parameters['phase'] * 2 ) % (2*np.pi) * 180 / np.pi
dist = injection_parameters['luminosity_distance']
GMST = bc_utils.gps_time_to_gmst(injection_parameters['geocent_time']) # is GMST !=0 possible in Yann's code??...
Rmin = 6.
m1min = 1.
m2min = 1.
f_low = minimum_frequency
sigpar, n = buf.CBC_Parameter_Allocation(m1, m2, Dec, RA, inc, psi, phic, dist, GMST, Rmin, m1min+m2min, f_low)
return sigpar, n
def injection_parameters_to_binaryFile(injection_parameters, minimum_frequency, binaryFileName):
m1 = injection_parameters['mass_1']
m2 = injection_parameters['mass_2']
Dec = injection_parameters['dec'] * 180 / np.pi
RA = injection_parameters['ra'] * 180 / np.pi
if 'iota' in injection_parameters.keys():
inc = injection_parameters['iota'] * 180 / np.pi
else:
inc = injection_parameters['theta_jn'] * 180 / np.pi
psi = injection_parameters['psi'] * 180 / np.pi
phic = injection_parameters['phase'] * 2 * 180 / np.pi
dist = injection_parameters['luminosity_distance']
GMST = bc_utils.gps_time_to_gmst(injection_parameters['geocent_time']) # is GMST !=0 possible in Yann's code??...
Rmin = 6.
m1min = 1.
m2min = 1.
f_low = minimum_frequency
m1max = 2.6
m2max = 2.6
Dmin = 1e-6
Dmax = 70
params = (RA, Dec, inc, psi, phic, dist, m1, m2, GMST, f_low, Rmin, m1min, m1max, m2min, m2max, Dmin, Dmax)
np.savetxt(binaryFileName, params)
return 0
def sigpar_to_ifo_parameters(sigpar):
"""
"""
duration = sigpar[16] # = Tobs = 1.1 * tchirp for us, cf buf.CBC_Parameter_Allocation()
sampling_frequency = 2.0 * sigpar[15]
# I don't use the + 2 here because it's kind of already done using Tobs = 1.1 * tchirp. Note sigpar[8] = tchirp = sigpar[16] / 1.1
# In bilby's example, the start_time is computed using the formula below where 'geocent_time' is the time of coalescence, ie our sigpar[8], but expressed in gps time when in the hmc sigpar[8] is counted from start_time
# start_time = injection_parameters['geocent_time'] + 2 - duration
# So we derive the start_time from our input GMST angle = sipar[18]. This will be the reference start_time through the whole algorithm, the geocent_time however will vary when sigpar[8] varies.
# However sigpar[18] is the GMST angle of the earth when the waveform arrives in Handford, what bilby wants is the start_time = time when the GW arrives at the center of the earth. Hence I should input my function `gmst_to_gps_time()` not with sigpar[18] but the GMST when the GW arrived at the earth center.
# bilby.gw.utils.time_delay_geocentric(H.vertex, np.array([0,0,0]), injection_parameters['ra'], injection_parameters['dec'], start_time)
H_time_delay_from_geocenter = 0.012139641777406887
omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
gmst_earth_center = sigpar[18] - H_time_delay_from_geocenter * omega_earth
start_time = gmst_to_gps_time(sigpar[18])
# start_time = geocent_time - sigpar[8]
# start_time = 1541967481.3621483#geocent_time - sigpar[8]
return duration, sampling_frequency, start_time
def gmst_to_gps_time(gmst):
"""
Converts a gmst angle to the corresponding gps_time closest to that of now.
"""
gps_time_now = time.time()
gmst_now = bc_utils.gps_time_to_gmst(gps_time_now)
gmst_diff = gmst_now - gmst
omega_earth = 2 * np.pi * (1 / 365.2425 + 1) / 86400.
number_of_seconds = gmst_diff / omega_earth
gps_time = gps_time_now - number_of_seconds
return gps_time
def sigpar_to_freq_and_indices(sigpar,n_freq):
"""
When computing waveforms and/or their derivatives, Fourier Transform purposes imposes that the array of frequencies has to start from zero and have a power of 2 length=n_freq. Hence the end frequency is higher than f_high.
We then set to zero all the waveform values corresponding to frequencies f not respecting : f_low <= f <= f_high; where f_low is an input of the binary.in file and f_high calculated from it.
This function returns:
return freqs, i_low, i_high
Where:
- freqs [array of floats]:
array of all the frequencies from 0 to n_freq*(delta_f-1) included
- i_low [int]:
s.t. freqs[i_low] is the 1st frequency where the waveform is non-zero
- i_high [int]:
s.t. freqs[i_high] is the last frequency where the waveform is non-zero
"""
delta_f = sigpar[17]
f_low = sigpar[13]
f_high = sigpar[14]
freqs = np.arange(0,n_freq*delta_f,delta_f)
i_low=int(f_low/delta_f)+1
i_high=int(f_high/delta_f)
return freqs, i_low, i_high
def sigpar_to_freq_and_indices_2(sigpar,n_freq):
"""
When computing waveforms and/or their derivatives, Fourier Transform purposes imposes that the array of frequencies has to start from zero and have a power of 2 length=n_freq. Hence the end frequency is higher than f_high.
We then set to zero all the waveform values corresponding to frequencies f not respecting : f_low <= f <= f_high; where f_low is an input of the binary.in file and f_high calculated from it.
This function returns:
return freqs, i_low, i_high
Where:
- freqs [array of floats]:
array of all the frequencies from 0 to n_freq*(delta_f-1) included
- i_low [int]:
s.t. freqs[i_low] is the 1st frequency where the waveform is non-zero
- i_high [int]:
s.t. freqs[i_high] is the last frequency where the waveform is non-zero
"""
duration, sampling_frequency, start_time = sigpar_to_ifo_parameters(sigpar)
import Library.set_injection_parameters as set_inj
duration = set_inj.ini_file_to_dict()['meta']['duration']
sampling_frequency = set_inj.ini_file_to_dict()['meta']['sampling_frequency']
freqs = bc_utils.create_frequency_series(sampling_frequency, duration)
f_low = sigpar[13]
f_high = sigpar[14]
delta_f = freqs[1] - freqs[0]
i_low=int(f_low/delta_f)+1
i_high=int(f_high/delta_f)
return freqs, i_low, i_high
def sigpar_to_n_freq_as_computed_in_bilby(sigpar):
"""
This comes from the computation in bilby.core.utils.create_frequency_series()
"""
duration, sampling_frequency, start_time = sigpar_to_ifo_parameters(sigpar)
number_of_samples = duration * sampling_frequency
number_of_samples = int(np.round(number_of_samples))
# prepare for FFT
number_of_frequencies = (number_of_samples - 1) // 2
if number_of_frequencies % 2 == 1:
number_of_frequencies += 2
else:
# no Nyquist frequency when N=odd
number_of_frequencies += 1
return number_of_frequencies
def get_next_random_numbers(n):
# global random_numbers
next_random_numbers = conf.random_numbers[0:n]
conf.random_numbers = conf.random_numbers[n:]
if n==1: next_random_numbers = next_random_numbers[0]
return next_random_numbers
def overlap_LVC_RiemannSum(CA1, CA2, shf, Tobs, n, flow, fhigh):
d1=buf.LVC_RiemannSum(CA1, CA1, shf, Tobs, n, flow, fhigh)
d2=buf.LVC_RiemannSum(CA2, CA2, shf, Tobs, n, flow, fhigh)
num=buf.LVC_RiemannSum(CA1, CA2, shf, Tobs, n, flow, fhigh)
return num/np.sqrt(d1*d2)
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