Skip to content

Uddeepta Deka / Updated Charged Lens

Go to a project
Commit c7c14fdd by Uddeepta Deka
Options

corrected lens equation

parent 859f04d7
Expand all
Inline Side-by-side
Showing with 38 additions and 82 deletions
notebooks/lens_class.ipynb
This source diff could not be displayed because it is too large. You can view the blob instead.
utils/lenses.py
......@@ -7,31 +7,25 @@ class ChargedLens(object):
Initialization params:
-------------------------
qeff : float
effective charge parameter
y1, y2 : float, float
dimensionless impact parameter
gtd_limit: float
max time delay gradient at images [default:4.5]
qeff : effective charge parameter
y1, y2 : dimensionless impact parameter
"""
def __init__(self, qeff, y1, y2=0., gtd_limit=4.5):
def __init__(self, qeff, y1, y2=0.):
"""
params:
---------------
qeff : float
qeff: float
effective charge
y1 : float
impact parameter
y2 : float [default:0.0]
impact parameter
gtd_limit : float
max time delay gradient at images [default:4.5]
"""
# TODO set qeff according to the values of Dl, Ds, Dls and Q
self.qeff = qeff
self.y1 = y1
self.y2 = y2
self.gtd_limit = gtd_limit
self.img, _ = self.imageLocations()
self.img, _ = self.imageLocations(mag_cut=True, mag_tol=1e-6)
self.img_midx = []
self.img_mag = []
self.img_td = []
......@@ -65,36 +59,11 @@ class ChargedLens(object):
def gradTimeDelay(self, x1, x2=0):
""" Gradient of time delay function """
a1, a2 = self.alpha(x1, x2)
dtdx1 = x1 - self.y1 - a1
dtdx2 = x2 - self.y2 - a2
x = np.sqrt(x1**2 + x2**2)
dtdx1 = x1 - self.y1 - (x1/x**2) + (self.qeff*x1/x**3)
dtdx2 = x2 - self.y2 - (x2/x**2) + (self.qeff*x2/x**3)
return dtdx1, dtdx2
def critLines(self):
""" Tangential critical points are mapped to y=0.
Radial critical points are obtained here"""
common_term = np.complex128(9 * self.qeff +
np.sqrt(3 * (1 + 27 * self.qeff**2)))
common_term = np.power(common_term, (1./3.))
expr1 = -1./(3**(1/3) * common_term) + (common_term/3**(2/3))
expr2 = (1 + 1j*np.sqrt(3))/(2*3**(1/3)*common_term) - (1 - 1j*np.sqrt(3))*common_term/(2*3**(2/3))
expr3 = (1 - 1j*np.sqrt(3))/(2*3**(1/3)*common_term) - (1 + 1j*np.sqrt(3))*common_term/(2*3**(2/3))
crit_pts = []
for crit_pt in [expr1, expr2, expr3]:
if (abs(np.imag(crit_pt)) < 1e-4):
crit_pts.append(np.real(crit_pt))
crit_pts = np.asarray(crit_pts)
return crit_pts
def caustic(self):
""" Radial caustic points """
crit_pts = self.critLines()
caustic_pts = []
for crit_pt in crit_pts:
caustic_pts.append(self.lensEquation(crit_pt))
caustic_pts = np.asarray(caustic_pts)
return caustic_pts
def detTrace(self, x1, x2):
""" Returns the determinant and trace of the time delay Hessian matrix """
x = np.sqrt(x1**2 + x2**2)
......@@ -104,7 +73,7 @@ class ChargedLens(object):
tr = 2 - self.qeff / x**3
return det, tr
def magnification(self, x1, x2):
def magnification(self, x1, x2=0):
""" Inverse determinant of the Hessian matrix """
det, tr = self.detTrace(x1, x2)
return 1. / det
......@@ -125,51 +94,38 @@ class ChargedLens(object):
return -1 # undefined
def lensEquation(self, x):
""" right hand side of the lens equation """
return x - (1/x) + (self.qeff/x**2)
return x - (1/x) + (self.qeff * x/(x**2)**1.5)
def critLines(self):
""" Tangential critical lines are mapped to y=0.
Radial critical lines are obtained here"""
common_term = np.complex128(9 * self.qeff +
np.sqrt(3 * (1 + 27 * self.qeff**2)))
common_term = np.power(common_term, (1./3.))
return -1./(3**(1/3) * common_term) + (common_term/3**(2/3))
def caustic(self):
x = self.critLines()
caustic_pts = x - 1/x + self.qeff/x**2
return caustic_pts
def imageLocations(self):
def imageLocations(self, gtd_cut=False, mag_cut=False, gtd_tol=1e-1, mag_tol=1e-3):
"""
Returns the location of images in the format
[[im1_x1, im2_x1, ...], [im1_x2, im2_x2, ...]]
"""
# for now we are doing it for y2 = 0
def eq1(x):
return x*x*self.y1
def eq2(x):
return x**3 - x + self.qeff
x = np.arange(-4, 4, 0.0005)
f1 = eq1(x)
f2 = eq2(x)
possible_img_locs = x[np.argwhere(np.diff(np.sign(f1 - f2))).flatten()]
x = np.arange(-2, 4, 0.0005)
possible_img_locs = x[np.argwhere(np.diff(np.sign(self.lensEquation(x) - self.y1))).flatten()]
## putting a cut on the grad time delay
if gtd_cut:
gtd1, gtd2 = self.gradTimeDelay(possible_img_locs, np.zeros_like(possible_img_locs))
gtdnorm = np.sqrt(gtd1**2 + gtd2**2)
idx = np.where(gtdnorm<self.gtd_limit)[0]
img_locs = np.array(possible_img_locs[idx])
return [img_locs, np.zeros_like(img_locs)]
def imageLocationsAnalytic(self, grad_td_check=True):
""" Analytically obtained locations of images """
val = 81 * self.qeff**2 - 3 * (4 + self.y1**2) - 6 * self.qeff * self.y1 * (9 + 2 * self.y1**2)
if np.sign(val)==-1:
common_term = np.complex128(27 * self.qeff - 9 * self.y1 - 2 * self.y1**3 + 3j * np.sqrt(np.abs(val)))
else:
common_term = np.complex128(27 * self.qeff - 9 * self.y1 - 2 * self.y1**3 + 3 * np.sqrt(val))
common_term = np.power(common_term, (1./3.))
expr1 = (1/6) * (2 * self.y1 - ((2**(4/3) * (3 + self.y1**2)) / common_term) - 2**(2/3) * common_term)
expr2 = (1/12) * (4 * self.y1 + ((2**(4/3) * (1 + 1j * np.sqrt(3)) * (3 + self.y1**2)) /
common_term) + 2**(2/3) * (1 - 1j * np.sqrt(3)) * common_term)
expr3 = (1/12) * (4 * self.y1 + ((2**(4/3) * (1 - 1j * np.sqrt(3)) * (3 + self.y1**2)) /
common_term) + 2**(2/3) * (1 + 1j * np.sqrt(3)) * common_term)
img_locs = []
for imgloc in [expr1, expr2, expr3]:
if (abs(np.imag(imgloc)) < 1e-4):
if grad_td_check:
gtd1, gtd2 = self.gradTimeDelay(imgloc, 0)
gtdnorm = np.sqrt(gtd1**2 + gtd2**2)
if gtdnorm < self.gtd_limit:
img_locs.append(np.real(imgloc))
else:
img_locs.append(np.real(imgloc))
img_locs = np.asarray(img_locs)
return img_locs
\ No newline at end of file
gtd_pass = gtdnorm < gtd_tol
possible_img_locs = possible_img_locs[gtd_pass]
## putting a magnification cutoff
if mag_cut:
img_mags = self.magnification(possible_img_locs, 0)
mag_pass = np.abs(img_mags) > mag_tol
possible_img_locs = possible_img_locs[mag_pass]
return [possible_img_locs, np.zeros_like(possible_img_locs)]
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment