Gwfish janosch

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,37 +394,26 @@ 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:
             i_final = np.searchsorted(frequencyvector, max_frequency_cutoff)
 
         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)

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

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()