introducing IMRPhenomD

Codes/set_injection_parameters.py

@@ -78,8 +78,20 @@ def ini_file_to_dict(file_path='../examples/GW170817.ini'):
     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)
+    approximant = 'TaylorF2'
+    if approximant == 'TaylorF2':
+        f_high = Compute_fmax(R_isco, m1, m2)
+        f_high_run = Compute_fmax(R_isco, m1_min, m2_min)
+    elif approximant = 'IMRPhenomD':
+        # cf:
+        # - https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom__c.html#gad3f98acfe9527259a7f73a8ef69a2f7b
+        # - line 103 of https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_i_m_r_phenom_d_8c_source.html#l00103
+        LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
+        f_CUT = 0.2
+        M_sec = (m1 + m2) * LAL_MTSUN_SI
+        M_min_sec = (m1_min + m2_min) * LAL_MTSUN_SI
+        f_high = f_CUT / M_sec
+        f_high_run = f_CUT / M_min_sec
     # import IPython; IPython.embed();
     duration = tc_3p5PN + 2
     # sampling_frequency = 4096

Library/bilby_waveform.py

@@ -13,10 +13,16 @@ import Library.config as conf
 def get_source_frame_polarizations(parameters, minimum_frequency, interferometers):
 
     # Fixed arguments passed into the source model. The analysis starts at 40 Hz.
-    approximant = 'TaylorF2'
+    # approximant = 'TaylorF2'
     # approximant = 'SpinTaylorF2'
     # approximant = 'IMRPhenomPv2'
-    waveform_arguments = dict(waveform_approximant=approximant, reference_frequency=0, minimum_frequency=minimum_frequency, maximum_frequency=0)
+    approximant = 'IMRPhenomD'
+    if approximant == 'TaylorF2':
+        maximum_frequency = 0
+    else:
+        maximum_frequency = 4096
+    # waveform_arguments = dict(waveform_approximant=approximant, reference_frequency=50, minimum_frequency=minimum_frequency)
+    waveform_arguments = dict(waveform_approximant=approximant, reference_frequency=0, minimum_frequency=minimum_frequency, maximum_frequency=maximum_frequency)
 
     # Create the waveform_generator using a LAL Binary Neutron Star source function
     waveform_generator = bilby.gw.WaveformGenerator(

Tests/steppy_hplus_over_m1_different_approx.py

+import numpy as np
+import matplotlib.pyplot as plt
+try:
+    import lal
+    import lalsimulation as lalsim
+except ImportError:
+    logger.warning("Please install lal and lalsimulation to run this script.")
+
+
+parsec = 3.0856775814671916e+16  # m
+solar_mass = 1.9884754153381438e+30  # Kg
+SOL = 2.99792458e8                                       # Speed of light / m/s
+G = 6.67428e-11
+GEOM = (G * solar_mass / (SOL*SOL*SOL))
+
+LAL_PI = 3.141592653589793238462643383279502884
+LAL_MSUN_SI = 1.988546954961461467461011951140572744e30
+LAL_MTSUN_SI = 4.925491025543575903411922162094833998e-6
+
+def Compute_fISCO(Rmin, m1, m2):
+    """
+    Compute the frequency at the Innermost Stable Circular Orbit
+    """
+
+    vISCO = np.sqrt(1.0 / Rmin)
+    Mt = m1 + m2
+    fmax = vISCO**3 / (np.pi * GEOM * Mt)      # fmax signal
+    return fmax
+
+def Compute_LAL_fISCO(m1_SI, m2_SI):
+    m1 = m1_SI / LAL_MSUN_SI
+    m2 = m2_SI / LAL_MSUN_SI
+    m_sec = (m1 + m2) * LAL_MTSUN_SI
+    piM = LAL_PI * m_sec
+    vISCO = 1 / np.sqrt(6)
+    fISCO = vISCO**3 / piM
+
+    return fISCO
+
+if __name__=='__main__':
+
+    luminosity_distance = 40
+    mass_1 = 1.47 + 160*1e-7
+    mass_2 = 1.27
+    spin_1x = 0.0
+    spin_1y = 0.0
+    spin_1z = 0.0
+    spin_2x = 0.0
+    spin_2y = 0.0
+    spin_2z = 0.0
+    iota = 2.81
+    phase = 5.497787143782138
+    longitude_ascending_nodes = 0.0
+    eccentricity = 0.0
+    mean_per_ano = 0.0
+    # delta_frequency = 0.017
+    delta_frequency = 0.016973449670253923
+    minimum_frequency = 30.0
+    reference_frequency = 0.0
+    # approximant_str = 'TaylorF2'
+    # approximant_str = 'IMRPhenomPv2'
+    approximant_str = 'IMRPhenomD'
+    approximant = lalsim.GetApproximantFromString(approximant_str)
+    # approximant = 5 # Corresponds to `TaylorF2`
+    if approximant_str == 'TaylorF2':
+        maximum_frequency = 0
+    else:
+        # maximum_frequency = 0
+        maximum_frequency = 4096
+
+    lambda_1 = 0.0
+    lambda_2 = 0.0
+
+
+    pn_amplitude_order = 0
+    pn_spin_order = -1
+    pn_tidal_order = -1
+    pn_phase_order = -1
+
+    if pn_amplitude_order != 0:
+        start_frequency = lalsim.SimInspiralfLow2fStart(
+            minimum_frequency, int(pn_amplitude_order), approximant)
+    else:
+        start_frequency = minimum_frequency
+
+    # waveform_dictionary = lal.CreateDict()
+    # lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(waveform_dictionary, int(pn_spin_order))
+    # lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(waveform_dictionary, int(pn_tidal_order))
+    # lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(waveform_dictionary, int(pn_phase_order))
+    # lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(waveform_dictionary, int(pn_amplitude_order))
+    # lalsim.SimInspiralWaveformParamsInsertTidalLambda1(waveform_dictionary, lambda_1)
+    # lalsim.SimInspiralWaveformParamsInsertTidalLambda2(waveform_dictionary, lambda_2)
+
+
+    luminosity_distance = luminosity_distance * 1e6 * parsec
+    mass_2 = mass_2 * solar_mass
+
+    step = 1e-7
+    nb_points = 8
+    mass_1_min = mass_1 - int(nb_points/2) * step
+    mass_1_max = mass_1 + int(nb_points/2) *step
+    masses_1 = np.linspace(mass_1_min, mass_1_max, nb_points)
+
+    steppy_sub_product = []
+    hplus_list = []
+    hcross_list = []
+    f_isco_LAL_list = []
+
+    for m1 in masses_1:
+
+        waveform_dictionary = lal.CreateDict()
+        lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(waveform_dictionary, int(pn_spin_order))
+        lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(waveform_dictionary, int(pn_tidal_order))
+        lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(waveform_dictionary, int(pn_phase_order))
+        lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(waveform_dictionary, int(pn_amplitude_order))
+        lalsim.SimInspiralWaveformParamsInsertTidalLambda1(waveform_dictionary, lambda_1)
+        lalsim.SimInspiralWaveformParamsInsertTidalLambda2(waveform_dictionary, lambda_2)
+
+        mass_1 = m1 * solar_mass
+
+        hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
+        mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, luminosity_distance, iota, phase,
+        longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
+        start_frequency, maximum_frequency, reference_frequency,
+        waveform_dictionary, approximant)
+
+        # breakpoint()
+
+        steppy_sub_product.append(hplus.data.data.real.sum())
+        hplus_list.append(hplus.data.data)
+        hcross_list.append(hcross.data.data)
+
+        # f_isco = Compute_fISCO(6, m1/solar_mass, mass_2/solar_mass)
+        # f_isco_LAL = Compute_LAL_fISCO(m1, mass_2)
+        # f_isco_LAL_list.append(f_isco_LAL)
+        # print('f_isco     = {}'.format(f_isco))
+        # print('f_isco_LAL = {} - array should be of size {} and be non-zero from iStart = {} \n'.format(f_isco_LAL, int(f_isco_LAL/delta_frequency + 1), int(minimum_frequency/delta_frequency)+1))
+        # Note: the same steppy behavior happenss with `hcross.data.data.imag.sum()` but does not with `hplus.data.data.imag.sum()` or `hcross.data.data.real.sum()`
+
+    steppy_sub_product_np = np.asarray(steppy_sub_product)
+    dhplus_real_sum = (steppy_sub_product_np[1:] - steppy_sub_product_np[:-1]) / step
+    # breakpoint()
+    # import IPython; IPython.embed()
+    # PLOT
+    from matplotlib import rc
+
+    plt.rcParams['text.usetex'] = False
+
+    plt.rcParams['figure.titlesize'] = 18
+    plt.rcParams['figure.titleweight'] = 'bold'
+
+    plt.rcParams['axes.titlesize'] = 15
+    plt.rcParams['axes.labelsize'] = 12
+    plt.rcParams['axes.titleweight'] = 'normal'
+
+    plt.rcParams['xtick.labelsize'] = 12
+    plt.rcParams['ytick.labelsize'] = 12
+    plt.rcParams['font.size'] = 10
+
+
+    fig0 = plt.figure(figsize=(16,8))
+    fig0.suptitle('Steppy behavior of source frame waveform when varying m1 - {}'.format(approximant_str))
+    # fig, axs = plt.subplots(2, 1, figsize=(16,8), gridspec_kw = {'height_ratios': [10, 1]})
+    ax0 = plt.subplot(111)
+    ax0.plot(masses_1, steppy_sub_product, linewidth=1, marker='o', markersize=0.5, label='hplus.real.sum()')
+    # ax0.plot(masses_1, steppy_sub_product, linewidth=1, marker='o', markersize=3, label='hplus.real.sum()')
+    ax0.set_xlabel('mass_1 (in solar masses)')
+    ax0.set_ylabel('hplus.real.sum()')
+    ax0.set_title('Step size between points =  {:.0e}'.format(step))
+    ax0.legend(loc=(0.85,0.7))
+
+    ax1 = ax0.twinx()
+    ax1.plot(masses_1[:-1], dhplus_real_sum, linewidth=1, marker='o', markersize=0.5, color='red', label='dhplus.real.sum()')
+    # ax1.plot(masses_1[:-1], dhplus_real_sum, linewidth=1, marker='o', markersize=3, color='red', label='dhplus.real.sum()')
+    ax1.tick_params('y', colors='red')
+    ax1.legend(loc=(0.85,0.65))
+
+
+    fig0.tight_layout()
+    # fig.text(x=0.4, y=0.93, s='Step size between points =  {:.0e}'.format(step))
+    fig0.subplots_adjust(top=0.88) # Needed so that fig.title and ax.title don't overlap
+
+    # file_name = './.steppy_hplus_over_m1/hplus_over_m1_{}.png'.format(approximant_str)
+    # fig0.savefig(file_name)
+
+
+
+
+
+    # nb_of_freq = len(hplus_list[0])
+    # frequency_array = np.linspace(0, nb_of_freq*delta_frequency, nb_of_freq)
+    #
+    # hplus_list[2] = np.append(hplus_list[2], 0)
+    # hplus_list[3] = np.append(hplus_list[3], 0)
+    #
+    # colors = ['blue', 'orange', 'green', 'red']
+    #
+    #
+    # fig1 = plt.figure(figsize=(16,8))
+    # fig1.suptitle('Steppy behavior of source frame waveform when varying one component mass')
+    # # fig, axs = plt.subplots(2, 1, figsize=(16,8), gridspec_kw = {'height_ratios': [10, 1]})
+    # ax1 = plt.subplot(111)
+    # for i, hplus in enumerate(hplus_list):
+    #     ax1.plot(frequency_array, hplus.real, linewidth=1, color=colors[i], marker='o', label='hplus.real - point #{}'.format(i+1))
+    #     ax1.axvline(x=f_isco_LAL_list[i], linewidth=0.5, color=colors[i], label='f_ISCO_LAL- point #{}'.format(i+1))
+    # ax1.set_xlabel('frequency')
+    # ax1.set_ylabel('hplus.real')
+    # ax1.set_title('Real part of hplus for different mass_1')
+    # ax1.legend(loc=(0.85,0.7))
+    #
+    # fig1.tight_layout()
+    # # fig.text(x=0.4, y=0.93, s='Step size between points =  {:.0e}'.format(step))
+    # fig1.subplots_adjust(top=0.88) # Needed so that fig.title and ax.title don't overlap
+    #
+    # # file_name = './.steppy_hplus_over_m1/steppy_hplus_over_m1.png'
+    # # fig.savefig(file_name)
+
+
+    plt.show()