code: def derivatives(self, x, y, theta_E, center_x=0, center_y=0):
          
          x_shift = x - center_x
          y_shift = y - center_y
          R = np.sqrt(x_shift*x_shift + y_shift*y_shift)
          if isinstance(R, int) or isinstance(R, float):
              a = theta_E / max(0.000001, R)
          else:
              a=np.empty_like(R)
              r = R[R > 0]  #in the SIS regime
              a[R == 0] = 0
              a[R > 0] = theta_E / r
          f_x = a * x_shift
          f_y = a * y_shift
          return f_x, f_y
docs: returns df/dx and df/dy of the function