& = & \frac{\theta_{E}}{2}\frac{\lambda(e)}{\sqrt{(1-e)^{-1}x^2+(1-e)y^2}},\\
\theta_{E} & = & 4\pi \Big(\frac{\sigma}{c}\Big)^{2}\frac{D_{\rm ls}}{D_{\rm s}},\end{aligned}$$ where $\theta_{E}$ is the Einstein radius, $\sigma$ is the line-of-sight velocity dispersion of the lens galaxy, $D_{\rm s}$ and $D_{\rm ls}$ are angular diameter distances, respectively, between the observer and the source, and the lens galaxy and the source. The convergence depends on the shape of the lens galaxy through the ellipticity $e$ of the projected lensing mass surface density, which includes the contributions from dark matter and baryonic components, both relevant for the statistics of strong lensing images [@Cas2018]. We assume that the ellipticity is given by the light distribution. Based on observational constraints, we adopt a Gaussian distribution for the ellipticity with a mean of 0.3 and dispersion of 0.16, with a truncation at 0.1 and 0.9 [@Ogu2008]. The function $\lambda(e)$ is the so-called dynamical normalization and it depends on the deprojected shape of the lens galaxy. Following @Chae2003, we assume that both oblate and prolate ellipsoids approximate the actual shapes of lens galaxies with equal probabilities, implying $\lambda(e)\approx1$.

In order to make the simulated lensing more realistic, we account for the effect of the lens environment by including external shear [@Koch1991; @Witt1997; @Kee1997] with a 