diff --git a/GWFish/modules/detection.py b/GWFish/modules/detection.py index d4aca87a..afb5a1c1 100644 --- a/GWFish/modules/detection.py +++ b/GWFish/modules/detection.py @@ -15,7 +15,7 @@ PSD_PATH = Path(__file__).parent.parent / 'detector_psd' # used when redefining the time and frequency vectors -N_FREQUENCY_POINTS = 1000 +N_FREQUENCY_POINTS = 2000 class DetectorComponent: @@ -143,9 +143,6 @@ def __init__(self, name: str, config=DEFAULT_CONFIG): fmax = eval(str(detector_def['fmax'])) spacing = str(detector_def['spacing']) - - - if spacing == 'linear': df = eval(str(detector_def['df'])) self.frequencyvector = np.linspace(fmin, fmax, int((fmax - fmin) / df) + 1) @@ -324,55 +321,58 @@ def projection(parameters, detector, polarizations, timevector, redefine_tf_vect f_max = parameters.get('max_frequency_cutoff', None) detector_lifetime = getattr(detector, 'mission_lifetime', None) - in_band_slice, new_timevector = in_band_window( + in_band_slice = in_band_window( np.squeeze(timevector), np.squeeze(detector.frequencyvector), detector_lifetime, - f_max, - redefine_timevector=redefine_tf_vectors + f_max ) - if redefine_tf_vectors: - new_fmin = detector.frequencyvector[in_band_slice.start - 1, 0] - new_fmax = detector.frequencyvector[in_band_slice.stop + 1, 0] - new_frequencyvector = np.geomspace(new_fmin, new_fmax, num=N_FREQUENCY_POINTS)[:, None] - temp_timevector = t_of_f_PN(parameters, new_frequencyvector) - in_band_slice, new_timevector = in_band_window( - np.squeeze(temp_timevector), - np.squeeze(new_frequencyvector), - detector_lifetime, - f_max, - redefine_timevector=True, - final_time=parameters['geocent_time'] - ) - - proj = np.zeros_like(new_timevector) + proj = np.zeros_like(timevector) if is_null_slice(in_band_slice): return proj with warnings.catch_warnings(): warnings.simplefilter('ignore', AstropyWarning) if detector.location == 'earth': - proj = projection_earth(parameters, detector, polarizations, new_timevector, in_band_slice, long_wavelength_approx = long_wavelength_approx) + proj = projection_earth(parameters, detector, polarizations, timevector, in_band_slice, long_wavelength_approx = long_wavelength_approx) elif detector.location == 'moon': - proj = projection_moon(parameters, detector, polarizations, new_timevector, in_band_slice) + proj = projection_moon(parameters, detector, polarizations, timevector, in_band_slice) elif detector.location == 'solarorbit': - proj = projection_solarorbit(parameters, detector, polarizations, new_timevector, in_band_slice) + proj = projection_solarorbit(parameters, detector, polarizations, timevector, in_band_slice) else: print('Unknown detector location') exit(0) - if redefine_tf_vectors: - return proj, new_timevector, new_frequencyvector return proj +def new_tf_vectors(parameters, detector, timevector, time_reset=True): + f_max = parameters.get('max_frequency_cutoff', None) + detector_lifetime = getattr(detector, 'mission_lifetime', None) + + in_band_slice = in_band_window( + np.squeeze(timevector), + np.squeeze(detector.frequencyvector), + detector_lifetime, + f_max + ) + + new_fmin = detector.frequencyvector[np.maximum(in_band_slice.start-2,0), 0] + new_fmax = detector.frequencyvector[np.minimum(in_band_slice.stop+2, len(detector.frequencyvector)-1), 0] + new_frequencyvector = np.geomspace(new_fmin, new_fmax, num=N_FREQUENCY_POINTS)[:, None] + new_timevector = t_of_f_PN(parameters, new_frequencyvector) + if time_reset: + time_shift = new_timevector[0, 0] + new_timevector -= time_shift + parameters['geocent_time'] = parameters['geocent_time']-time_shift + + return new_timevector, new_frequencyvector, parameters def in_band_window( timevector, frequencyvector, detector_lifetime, max_frequency_cutoff, - redefine_timevector=False, final_time=None ): """Truncate the evolution of the source to the detector lifetime. @@ -394,12 +394,11 @@ def in_band_window( if max_frequency_cutoff is None: i_final = len(timevector) fmax_time = final_time - new_timevector = timevector else: if max_frequency_cutoff <= frequencyvector[0]: warnings.warn("The max_frequency given is lower than the lowest frequency for the detector." "Returning a zero projection.") - return slice(0, 0), timevector + return slice(0, 0) elif max_frequency_cutoff >= frequencyvector[-1]: i_final = -1 else: @@ -407,24 +406,14 @@ def in_band_window( fmax_time = timevector[i_final] - if redefine_timevector: - # potentially dangerous loss of precision here... - # shift the timevector so that the cutoff is placed at the given gps time - new_timevector = timevector + final_time - fmax_time - else: - new_timevector = timevector - if detector_lifetime is None: i_initial = 0 - elif final_time - detector_lifetime < new_timevector[0]: + elif final_time - detector_lifetime < timevector[0]: i_initial = 0 else: - if redefine_timevector: - i_initial = np.searchsorted(new_timevector, final_time - detector_lifetime) - else: - i_initial = np.searchsorted(new_timevector, fmax_time - detector_lifetime) + i_initial = np.searchsorted(timevector, fmax_time - detector_lifetime) - return slice(i_initial, i_final), new_timevector + return slice(i_initial, i_final) def projection_solarorbit(parameters, detector, polarizations, timevector, in_band_slice=slice(None)): ff = detector.frequencyvector[in_band_slice] @@ -477,7 +466,7 @@ def projection_earth(parameters, detector, polarizations, timevector, in_band_sl # timevector = parameters['geocent_time'] * np.ones_like(timevector) # switch off Earth's rotation - nf = len(polarizations[:, 0]) + nf = len(timevector) ff = detector.frequencyvector[in_band_slice, :] components = detector.components @@ -598,7 +587,7 @@ def projection_moon(parameters, detector, polarizations, timevector, in_band_sli # timevector = parameters['geocent_time'] * np.ones_like(timevector) # switch off Earth's rotation - nt = len(polarizations[:, 0]) + nt = len(timevector) components = detector.components proj = np.zeros((nt, len(components)), dtype=complex) diff --git a/GWFish/modules/horizon.py b/GWFish/modules/horizon.py index 09c91cc6..961309b5 100644 --- a/GWFish/modules/horizon.py +++ b/GWFish/modules/horizon.py @@ -15,17 +15,16 @@ from astropy.cosmology import Planck18 import astropy.cosmology as cosmology -import astropy.units as u -from scipy.optimize import brentq, minimize, dual_annealing +import scipy.optimize as opt -from .detection import SNR, Detector, projection, Network +from .detection import SNR, Detector, projection, Network, new_tf_vectors from .waveforms import LALFD_Waveform, DEFAULT_WAVEFORM_MODEL, Waveform DEFAULT_RNG = np.random.default_rng(seed=1) -MIN_REDSHIFT = 1e-20 -MAX_REDSHIFT = 1e6 +MIN_REDSHIFT = 1e-30 +MAX_REDSHIFT = 5000 def compute_SNR( params: "Union[dict[str, float], pd.DataFrame]", @@ -44,6 +43,8 @@ def compute_SNR( :return: the SNR """ + detector_frequency_vector_keeper = detector.frequencyvector + data_params = { 'frequencyvector': detector.frequencyvector, 'f_ref': 50. @@ -52,15 +53,21 @@ def compute_SNR( polarizations = waveform_obj() timevector = waveform_obj.t_of_f - args = (params, detector, polarizations, timevector) - if redefine_tf_vectors: - signal, timevector, frequencyvector = projection(*args, redefine_tf_vectors=True) - else: - signal = projection(*args) - frequencyvector = detector.frequencyvector + timevector, frequencyvector, params = new_tf_vectors(params, detector, timevector, time_reset=True) + detector.frequencyvector = frequencyvector + data_params['frequencyvector'] = frequencyvector + + waveform_obj = waveform_class(waveform_model, params, data_params) + polarizations = waveform_obj() + + args = (params, detector, polarizations, timevector) + signal = projection(*args) + + component_SNRs = SNR(detector, signal, frequencyvector=np.squeeze(detector.frequencyvector)) + + detector.frequencyvector = detector_frequency_vector_keeper - component_SNRs = SNR(detector, signal, frequencyvector=np.squeeze(frequencyvector)) return np.sqrt(np.sum(component_SNRs**2)) def compute_SNR_network( @@ -129,7 +136,7 @@ def SNR_error(redshift): warnings.warn('The source is completely out of band') return 0., 0. - redshift, r = brentq(SNR_error, MIN_REDSHIFT, MAX_REDSHIFT, full_output=True) + redshift, r = opt.brentq(SNR_error, MIN_REDSHIFT, MAX_REDSHIFT, full_output=True) if not r.converged: raise ValueError('Horizon computation did not converge!') @@ -209,11 +216,11 @@ def to_minimize(x): if 'maxiter' not in minimizer_kwargs: minimizer_kwargs['maxiter'] = 100 - res = dual_annealing( + res = opt.dual_annealing( func=to_minimize, bounds=[ (0, 2*np.pi), - (-np.pi, np.pi), + (-np.pi/2., np.pi/2.), ], x0=x0, **minimizer_kwargs diff --git a/GWFish/modules/waveforms.py b/GWFish/modules/waveforms.py index 77609eb2..69ee9364 100644 --- a/GWFish/modules/waveforms.py +++ b/GWFish/modules/waveforms.py @@ -992,7 +992,7 @@ def plot(self, output_folder='./'): plt.legend(fontsize = 8) plt.grid(which = 'both', color = 'lightgray', alpha = 0.5, linestyle = 'dashed', linewidth = 0.5) plt.xlabel('Frequency [Hz]') - plt.ylabel('$\phi$_prime') + plt.ylabel(r'$\phi$_prime') plt.savefig(output_folder + 'psi_prime_phenomD.png') plt.close() @@ -1012,7 +1012,7 @@ def plot(self, output_folder='./'): plt.legend(fontsize = 8) plt.grid(which = 'both', color = 'lightgray', alpha = 0.5, linestyle = 'dashed', linewidth = 0.5) plt.xlabel('Frequency [Hz]') - plt.ylabel('$\phi$') + plt.ylabel(r'$\phi$') plt.savefig(output_folder + 'psi_phenomD_zoomed.png') plt.close()