potential given by $$V(x,y)=\frac{\gamma}{2}(x^{2}-y^{2})\cos(2\theta_{\gamma})+\gamma xy\sin(2\theta_{\gamma}),$$ where $\gamma$ is the magnitude of the external shear and $\theta_{\gamma}$ is its position angle on the image plane. We assume that $\gamma$ follows a log-normal distribution with mean 0.05 and dispersion 0.2 dex, as expected for the external shear around early-type galaxies, calculated by ray tracing in cosmological simulations of the standard cosmological model [@Hol2003]. We expect that the employed simulation-based calibration of the external shear to a large extent accounts for lensing effects from realistic structures around lens elliptical galaxies, although the presence of a galaxy cluster may in addition affect the image configuration. The external shear in our model is uncorrelated with the orientation of the lens galaxy, i.e. it has a random orientation in the image plane.

We model the mass function of lens galaxies in terms of the velocity dispersion function of early-type galaxies, which are most common lens galaxies. The velocity dispersion function is well approximated by a modified Schechter function of the following form: $$\frac{dn}{d\sigma}=\phi_{\star}\Big(\frac{\sigma}{\sigma_{\star}}\Big)^{\alpha}
\exp\Big[-\Big(\frac{\sigma}{\sigma_{\star}}\Big)^{\beta}\Big]\frac{\beta}{\Gamma(\alpha/\beta)}\frac{1}{\sigma},$$ where $n$ is the comoving density of galaxies. We use the above probability density to draw random realizations of lens galaxies in our calculations. We adopt parameters derived from fitting this model 