updated with new hyp1f1 file

scripts/isolated_lens_ampfac.py

@@ -4,48 +4,23 @@ 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)
+def lens_ampl_function_wave_optics(w, y_l, intrp_hyp1f1=None):
+    """ frequency dependent lensing magnification """
+    x_m = (y_l + np.sqrt(y_l**2 +4)) / 2.
+    phi_m = (x_m - y_l)**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.vectorize(mpmath.hyp1f1)(0.5*w*1j, 1., 0.5*w*y_l**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))
+        hyp_fn = 10**intrp_hyp1f1['log_abs'](w, y_l)*np.exp(1j*intrp_hyp1f1['arg'](w, y_l))
     
     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
-    """
-            
+def lens_ampl_function_geom_optics(w, y): 
+    """ frequency dependent lensing magnification (geometric optics limit) """
     # 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 
@@ -57,30 +32,30 @@ def ampfac_geom(w, y):
     # 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 
+
+def lens_ampl_function(w, y_l, intrp_hyp1f1=None):
+    """ frequency dependent lensing magnification: hybrid using wave optics and geometric optics """
+
+    # if y > 3 geometric optics provides a good approximation 
+    if y_l > 3:
+        F_f = lens_ampl_function_geom_optics(w, y_l)
+
+    # if y <= 3, use a combination of wave optics (low freq) and geom optics (high freq)    
+    else:
+        w_max = 150
+        w_1 = w[np.where(w <= w_max)[0]]
+        w_2 = w[np.where(w > w_max)[0]]
+        
+        if len(w_1) == 0:
+            F_f = lens_ampl_function_geom_optics(w, y_l)
+        
+        else:
+            if len(w_2) != 0:
+                F_w = lens_ampl_function_wave_optics(w, y_l, intrp_hyp1f1=intrp_hyp1f1)
+                F_g = lens_ampl_function_geom_optics(w, y_l)
+                F_f = np.concatenate((F_w,F_g))
+
+            else:
+                F_f = lens_ampl_function_wave_optics(w, y_l, intrp_hyp1f1=intrp_hyp1f1)
+
+    return F_f