code: def fermat_potential(self, x_image, y_image, x_source, y_source, kwargs_lens, k=None):
          
  
          potential = self.potential(x_image, y_image, kwargs_lens, k=k)
          geometry = ((x_image - x_source)**2 + (y_image - y_source)**2) / 2.
          return geometry - potential
docs: fermat potential (negative sign means earlier arrival time)
  
          :param x_image: image position
          :param y_image: image position
          :param x_source: source position
          :param y_source: source position
          :param kwargs_lens: list of keyword arguments of lens model parameters matching the lens model classes
          :return: fermat potential in arcsec**2 without geometry term (second part of Eqn 1 in Suyu et al. 2013) as a list