[Python] GWFish's signal generation take much longer when expanding frequency range beyond...

I’m using the GWFish tutorial to simulate a binary neutron star (BNS) signal for the Einstein Telescope (ET). By default, the signal is projected onto the detector's sensitivity curve (defined in detectors.yaml with fmin=2 and fmax=2048 Hz) as the plot below: enter image description here

However, I need the signal before it enters ET’s detection band, so I tried two approaches:

1. 
   Custom frequency range: Used get_fd_signal_custom_range() to generate the signal over a broader range (f_min=0.0001, f_max=10000).
   
   
2. 
   Modified detectors.yaml: Manually expanded fmin/fmax in the ET configuration.
   

Both methods work but now take ~1 hour to run (vs. seconds for the default range). Additionally, the resulting characteristic strain shows oscillations instead of the expected quasi-straight line for a GW signal. Here's my script and the plot:

# In the file '/Users/.../myenv/lib/python3.9/site-packages/GWFish/detectors.yaml', for ET, I changed the values:
# fmin: 2.
# fmax: 2048.

# to:

# fmin: 0.0001
# fmax: 100000

# to calculate the signal in a broader range (out of the range of the sensitivity curve of ET)"


# suppress warning outputs for using lal in jupuyter notebook
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import GWFish.modules as gw
from tqdm import tqdm
import matplotlib
import matplotlib.pyplot as plt
import corner
import numpy as np
import pandas as pd
import json
import os
from astropy.cosmology import Planck18


# Event's parameters should be passed as Pandas dataframe

parameters = {
'chirp_mass': np.array([1.1858999987203738]) * (1 + 0.00980),
'mass_ratio': np.array([0.8308538032620448]),
'luminosity_distance': Planck18.luminosity_distance(0.00980).value,
'theta_jn': np.array([2.545065595974997]),
'ra': np.array([3.4461599999999994]),
'dec': np.array([-0.4080839999999999]),
'psi': np.array([0.]),
'phase': np.array([0.]),
'geocent_time': np.array([1187008882.4]),
'a_1':np.array([0.005136138323169717]),
'a_2':np.array([0.003235146993487445]),
'lambda_1':np.array([368.17802383555687]),
'lambda_2':np.array([586.5487031450857])}
parameters = pd.DataFrame(parameters)
parameters


# We choose a waveform approximant suitable for BNS analysis
# In this case we are taking into account tidal polarizability effects
waveform_model = 'IMRPhenomD_NRTidalv2'
f_ref = 10. # Reference frequency or 10. or 50.


# To access the available detectors, use the following command
gw.utilities.get_available_detectors()


# Choose the detector onto which you want to project the signal
detector = 'ET'

# The following function outputs the signal projected onto the chosen detector
signal, _ = gw.utilities.get_fd_signal(parameters, detector, waveform_model, f_ref) # waveform_model and f_ref are passed together

# I created a customized function in the module to include broader frequencies:
custom_signal, _ = gw.utilities.get_fd_signal_custom_range(parameters, detector, waveform_model, f_ref, f_min=0.0001, f_max=10000, f_bins=1000)


# Get the frequency vector from frequencyvector[:, 0] which is based on the detector's PSD.
frequency = gw.detection.Detector(detector).frequencyvector[:, 0]

#customize the frequency range:
custom_frequency = np.logspace(np.log10(1.e-4), np.log10(1.e4), num=1000) # Example range

# Assuming 'frequency' and 'signal' are already defined
characteristic_strain_1 = 2 * frequency * np.abs(signal[:, 0])
characteristic_strain_2 = 2 * custom_frequency * np.abs(custom_signal[:, 0])

# Apply a cut: Only plot where characteristic_strain < 1e-21
mask_1 = characteristic_strain_1 > 1e-23
mask_2 = characteristic_strain_2 > 1e-23

# add the detector's sensitivity curve and plot the characteristic strain
psd_data = gw.utilities.get_detector_psd(detector)
# psd_data is a dictionary with the keys 'frequency' and 'psd'

plt.loglog(custom_frequency[mask_2], characteristic_strain_2[mask_2], label='customed_%s' %waveform_model, color='pink')

#plt.loglog(frequency, 2 * frequency * np.abs(signal[:, 0]), label='%s' %waveform_model, color='red')

plt.loglog(frequency[mask_1], characteristic_strain_1[mask_1], label='%s' % waveform_model, color='red')

plt.loglog(psd_data[:, 0], np.sqrt(psd_data[:, 0] * psd_data[:, 1]), label='%s sensitivity' %detector, color='green')

plt.legend()
plt.xlabel('Frequency [Hz]')
plt.ylabel('Characteristic Strain')
plt.grid(linestyle='dotted', linewidth='0.6', which='both')
plt.title('Projected signal on %s' %detector)
plt.savefig('test.png')
plt.show()


enter image description here

Questions:

• 
  Why does expanding the frequency range slow down the computation so drastically? Is this due to the waveform model (IMRPhenomD_NRTidalv2) or PSD calculations? How to fix it?
  
  
• 
  Why does the characteristic strain oscillate at frequencies outside the detector’s sensitive band? How to fix it to have a classic straight line signal instead?
  

Continue reading...