set_injection_parameters.py 3.81 KB
Newer Older
Marc Arene's avatar
Marc Arene committed
1 2 3
import sys
# cf https://docs.python.org/3/tutorial/modules.html
sys.path.append('../')
4
import numpy as np
5 6
from configparser import ConfigParser

7
import bilby.gw.conversion as conv
8

9
from Headers.PN_Coefficients import * # Headers.Constants already imported in PN_Coeff
10
import Library.python_utils as pu
11
import Library.param_utils as paru
Marc Arene's avatar
Marc Arene committed
12

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56


def ComputeChirpTime3p5PN(f_low, m1, m2):
    """
    Computes chirp times at 3.5 PN.

    Inputs:
    -------
        f_low: [float] low frequency cutoff of the detector
        sigpar: [array_like] parameters of the system
    Outputs:
    --------

    """
    Mt = m1 + m2
    Mt_sec = Mt * GEOM           # Total mass in s
    etaM = m1 * m2 / Mt**2
    etaM2 = etaM*etaM

    Tau_0 = tf1 * Mt_sec / etaM
    Tau_2 = tf2 + tf3 * etaM
    Tau_3 = tf4 * PI
    Tau_4 = tf5 + tf6 * etaM + tf7 * etaM2
    Tau_5 = PI*(tf9 + tf8*etaM)
    Tau_6 = tf11 + tf12*PI2 + tf10*(EULER + np.log(4.0)) + ( tf13 + tf14*PI2 )*etaM + tf15*etaM2 + tf16*etaM2*etaM
    Tau_6_log = tf10
    Tau_7 = PI*(tf19 + tf18*etaM + tf17*etaM2 )
    v_0 = pow(PI*Mt_sec*f_low,1.0/3.0)
    v_0_2 = v_0*v_0
    v_0_3 = v_0_2*v_0
    v_0_4 = v_0_3*v_0
    v_0_5 = v_0_4*v_0
    v_0_6 = v_0_5*v_0
    v_0_7 = v_0_6*v_0

    Tau_chirp = Tau_0*pow(v_0,-8.0)*(1.0 + Tau_2*v_0_2 + Tau_3*v_0_3 + Tau_4 *v_0_4 +  Tau_5 *v_0_5 +  (Tau_6 + Tau_6_log*np.log(v_0))*v_0_6 +  Tau_7 *v_0_7 )
    return Tau_chirp

def Compute_fmax(Rmin, m1, m2):
    vmax = np.sqrt(1.0 / Rmin)
    Mt = m1 + m2
    fmax = vmax**3 / (PI * GEOM * Mt)      # fmax signal
    return fmax

57 58 59 60 61 62 63 64
def ini_file_to_dict(file_path='../examples/GW170817.ini'):

    config = ConfigParser()
    config.read(file_path)
    injection_parameters = pu.config_parser_to_dict(config)['parameters']


    injection_parameters = conv.generate_mass_parameters(injection_parameters)
65 66 67
    # injection_parameters['reduced_mass'] = (injection_parameters['mass_1'] * injection_parameters['mass_2']) / (injection_parameters['mass_1'] + injection_parameters['mass_2'])
    injection_parameters['reduced_mass'] = paru.component_masses_to_reduced_mass(injection_parameters['mass_1'], injection_parameters['mass_2'])

68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101


    # Set the duration and sampling frequency of the data segment that we're going
    # to inject the signal into. For the
    # TaylorF2 waveform, we cut the signal close to the isco frequency
    minimum_frequency=30.0
    m1 = injection_parameters['mass_1']
    m2 = injection_parameters['mass_2']
    tc_3p5PN = ComputeChirpTime3p5PN(minimum_frequency, m1, m2)
    # We integrate the signal up to the frequency of the "Innermost stable circular orbit (ISCO)"
    R_isco = 6.      # Orbital separation at ISCO, in geometric units. 6M for PN ISCO; 2.8M for EOB
    m1_min = 1
    m2_min = 1
    f_high = Compute_fmax(R_isco, m1, m2)
    f_high_run = Compute_fmax(R_isco, m1_min, m2_min)
    # import IPython; IPython.embed();
    duration = tc_3p5PN + 2
    # sampling_frequency = 4096
    # sampling_frequency = 4397
    sampling_frequency = 2 * f_high_run
    start_time = injection_parameters['geocent_time'] - tc_3p5PN

    meta_params = dict(
        minimum_frequency = minimum_frequency,
        tc_3p5PN = tc_3p5PN,
        f_high = f_high,
        f_high_run = f_high_run,
        duration = duration,
        sampling_frequency = sampling_frequency,
        start_time = start_time
    )

    injection_parameters['meta'] = meta_params

Marc Arene's avatar
Marc Arene committed
102 103

    # These two lines have been added to be able to plot the injected value on phi_c on the corner plots. Otherwise the `phi_c` and `tc_3p5PN` keys are not used in the code.
104
    injection_parameters['phi_c'] = (injection_parameters['phase'] * 2) % (2*np.pi)
Marc Arene's avatar
Marc Arene committed
105 106 107
    injection_parameters['tc_3p5PN'] = tc_3p5PN


108
    return injection_parameters
109 110 111 112 113 114 115 116 117




if __name__=='__main__':

    injection_parameters = ini_file_to_dict()

    pu.print_dict(injection_parameters, indent=3, align_keys=True)