fixed window_frac to 0.05

scripts/generate_ampfac_charged_lens.py

@@ -145,6 +145,7 @@ class ReconstructPeak:
     
     
     def find_peak_window(self):
     def find_peak_window(self):
         t_peak = self.Tj()
         t_peak = self.Tj()
+        assert t_peak > 1e-10
         t_start = (1- self.peak_window_frac) * t_peak
         t_start = (1- self.peak_window_frac) * t_peak
         t_end = (1 + self.peak_window_frac) * t_peak
         t_end = (1 + self.peak_window_frac) * t_peak
         start_of_window = np.where(self.t >= t_start)[0][0]
         start_of_window = np.where(self.t >= t_start)[0][0]
@@ -296,20 +297,6 @@ def main():
     w_max = 20
     w_max = 20
     dt = np.pi / w_max / dt_fac
     dt = np.pi / w_max / dt_fac
         
         
-    ## determining peak correction window
-    if y < 0.015:
-        win_frac = 0.5
-    elif y < 0.025:
-        win_frac = 0.4
-    elif y < 0.1:
-        win_frac = 0.2
-    elif y < 0.5:
-        win_frac = 0.05
-    elif y < 0.8:
-        win_frac = 0.02
-    elif y >=0.8:
-        win_frac = 0.01
-    
     tStart = time.perf_counter()
     tStart = time.perf_counter()
     print(f"Generating data for y = {y}, q={q} using {method} method ... ")
     print(f"Generating data for y = {y}, q={q} using {method} method ... ")
     # Creating the lens object
     # Creating the lens object
@@ -323,7 +310,7 @@ def main():
         t_dense = np.linspace(min(t_vals), max(t_vals), num)
         t_dense = np.linspace(min(t_vals), max(t_vals), num)
         Ft_interped = np.interp(t_dense, t_vals, Ft_vals)
         Ft_interped = np.interp(t_dense, t_vals, Ft_vals)
         # applying peak correction
         # applying peak correction
-        rp = ReconstructPeak(lens, t_dense, Ft_interped, peak_window_frac=win_frac)
+        rp = ReconstructPeak(lens, t_dense, Ft_interped, peak_window_frac=0.05)
         Ft_recons = rp.get_peak_reconstructed_Ft()
         Ft_recons = rp.get_peak_reconstructed_Ft()
         # doing the FFT
         # doing the FFT
         w_vals, Fw_vals = ampFac(t_dense, Ft_interped)
         w_vals, Fw_vals = ampFac(t_dense, Ft_interped)

scripts/isolated_lens_ampfac.py

+import numpy as np
+import time
+from scipy import special
+from mpmath import hyp1f1
+import pickle
+
+def ampfac_analytic(w, y, intrp_hyp1f1=None):
+    """ 
+    Computes the analytic amplification factor for isolated point mass lens
+    Input params:
+    ---------------
+    w : dimensionless angular frequency vector
+    y : impact parameter
+    
+    Returns:
+    ---------------
+    Complex amplification factor
+    """
+        
+    x_m = (y + np.sqrt(y**2 +4)) / 2.
+    phi_m = (x_m - y)**2/2. - np.log(x_m)
+    
+    if intrp_hyp1f1==None:
+        # evaluate the hypergeometric fun using mpmath package 
+        hyp_fn = np.vectorize(hyp1f1)(0.5*w*1j, 1., 0.5*w*y**2*1j, maxterms=1e6)
+        hyp_fn = np.array(hyp_fn.tolist(),dtype=complex)  # convert to numpy array from mpmath format 
+    else: 
+        # evaluate the interpolated function 
+        hyp_fn = 10**intrp_hyp1f1['log_abs'](w, y)*np.exp(1j*intrp_hyp1f1['arg'](w,y))
+    
+    return np.exp(np.pi*w/4.+0.5*w*1j*(np.log(w/2.)-2*phi_m))*special.gamma(np.ones(len(w))-0.5*w*1j)*hyp_fn
+
+def ampfac_geom(w, y): 
+    """ 
+    Amplification factor in geometric optics limit
+    for unresolved images
+    Eqn 15, Takahashi, Nakamura 2003
+    
+    Input params:
+    ---------------
+    w : dimensionless angular frequency vector
+    y : impact parameter
+
+    Returns:
+    ---------------
+    Complex amplification factor
+    """
+            
+    # time delay 
+    sqrt_ysqr_p4 = np.sqrt(y**2+4)
+    delta_Td_hat = 0.5*y*sqrt_ysqr_p4 + np.log((sqrt_ysqr_p4+y)/(sqrt_ysqr_p4-y)) # this is dimensionless 
+
+    # magnification 
+    mu_p = 0.5+(y**2+2)/(2*y*sqrt_ysqr_p4)
+    mu_m = 0.5-(y**2+2)/(2*y*sqrt_ysqr_p4)
+
+    # magnification function 
+    return np.sqrt(np.abs(mu_p)) - 1j*np.sqrt(np.abs(mu_m))*np.exp(1j*w*delta_Td_hat)
+
+def ampfac_hybrid(w, y, intrp_hyp1f1=None):
+    """ 
+    Hybrid amplification factor
+    
+    Input params:
+    ---------------
+    w : dimensionless angular frequency vector
+    y : impact parameter
+    
+    Returns:
+    ---------------
+    Complex amplification factor
+    """
+    
+    # if y >= 1.5 geometric optics provides a good approximation 
+    if y >= 2: 
+        F_f = ampfac_geom(w, y)
+    # if y < 2, use a combination of wave optics (low freq) and geom optics (high freq)
+    else: 
+        F_f=np.zeros(len(w),dtype='complex128')
+        wo_idx = w < 100
+        if len(w[wo_idx]) > 1: 
+            F_f[wo_idx] = ampfac_analytic(w[wo_idx], y, intrp_hyp1f1=intrp_hyp1f1)
+        if len(w[~wo_idx]) > 1:
+            F_f[~wo_idx] = ampfac_analytic(w[~wo_idx], y)
+            
+    return F_f 

utils/helper_functions.py

@@ -82,15 +82,16 @@ def lensed_fd_waveform(**params):
     lensed_hc = lensed_hc.cyclic_time_shift(3)
     lensed_hc = lensed_hc.cyclic_time_shift(3)
     return f, lensed_hp, lensed_hc
     return f, lensed_hp, lensed_hc
 
 
-def CalcLogMisMatch(xf1, xf2, Sh, f_min, f_max):
+def CalcLogMisMatch(xf1, xf2, Sh, f_min, f_max, optimized_match=False):
+    if optimized_match:
+        match, shift_idx = pycbc.filter.matchedfilter.optimized_match(
+            xf1, xf2, psd=Sh, low_frequency_cutoff=f_min, high_frequency_cutoff=f_max,
+            v1_norm=None, v2_norm=None,
+        )
+    else:
         match, shift_idx = pycbc.filter.matchedfilter.match(
         match, shift_idx = pycbc.filter.matchedfilter.match(
-        xf1,
-        xf2,
-        psd=Sh,
-        low_frequency_cutoff=f_min,
-        high_frequency_cutoff=f_max,
-        v1_norm=None,
-        v2_norm=None,
+            xf1, xf2, psd=Sh, low_frequency_cutoff=f_min, high_frequency_cutoff=f_max,
+            v1_norm=None, v2_norm=None,
         )
         )
     if match >= 1:
     if match >= 1:
         return -100
         return -100