diff --git a/pde/grids/base.py b/pde/grids/base.py
index a4a889a4..67be145e 100644
--- a/pde/grids/base.py
+++ b/pde/grids/base.py
@@ -512,55 +512,111 @@ def uniform_cell_volumes(self) -> bool:
         else:
             return all(np.asarray(vols).ndim == 0 for vols in self.cell_volume_data)
 
-    def difference_vector_real(self, p1: np.ndarray, p2: np.ndarray) -> np.ndarray:
+    def _difference_vector(
+        self,
+        p1: np.ndarray,
+        p2: np.ndarray,
+        *,
+        coords: CoordsType,
+        periodic: Sequence[bool],
+        axes_bounds: tuple[tuple[float, float], ...] | None,
+    ) -> np.ndarray:
         """return Cartesian vector(s) pointing from p1 to p2
 
         In case of periodic boundary conditions, the shortest vector is returned.
 
-        Warning:
-            This function assumes a Cartesian coordinate system and might thus not work
-            for curvilinear coordinates.
-
         Args:
             p1 (:class:`~numpy.ndarray`):
                 First point(s)
             p2 (:class:`~numpy.ndarray`):
                 Second point(s)
+            coords (str):
+                The coordinate system in which the points are specified.
+            periodic (sequence of bool):
+                Indicates which cartesian axes are periodic
+            axes_bounds (sequence of pair of floats):
+                Indicates the bounds of the cartesian axes
 
         Returns:
             :class:`~numpy.ndarray`: The difference vectors between the points with
-            periodic  boundary conditions applied.
+            periodic boundary conditions applied.
         """
-        diff = np.atleast_1d(p2) - np.atleast_1d(p1)
-        assert diff.shape[-1] == self.num_axes
+        x1 = self.transform(p1, source=coords, target="cartesian")
+        x2 = self.transform(p2, source=coords, target="cartesian")
+        if axes_bounds is None:
+            axes_bounds = self.axes_bounds
 
-        for i, periodic in enumerate(self.periodic):
-            if periodic:
-                size = self.axes_bounds[i][1] - self.axes_bounds[i][0]
+        diff = np.atleast_1d(x2) - np.atleast_1d(x1)
+        assert diff.shape[-1] == self.dim
+
+        for i, per in enumerate(periodic):
+            if per:
+                size = axes_bounds[i][1] - axes_bounds[i][0]
                 diff[..., i] = (diff[..., i] + size / 2) % size - size / 2
         return diff  # type: ignore
 
-    def distance_real(self, p1: np.ndarray, p2: np.ndarray) -> float:
+    def difference_vector(
+        self, p1: np.ndarray, p2: np.ndarray, *, coords: CoordsType = "grid"
+    ) -> np.ndarray:
+        """return Cartesian vector(s) pointing from p1 to p2
+
+        In case of periodic boundary conditions, the shortest vector is returned.
+
+        Args:
+            p1 (:class:`~numpy.ndarray`):
+                First point(s)
+            p2 (:class:`~numpy.ndarray`):
+                Second point(s)
+            coords (str):
+                The coordinate system in which the points are specified. Valid values are
+                `cartesian`, `cell`, and `grid`; see :meth:`~pde.grids.base.GridBase.transform`.
+
+        Returns:
+            :class:`~numpy.ndarray`: The difference vectors between the points with
+            periodic boundary conditions applied.
+        """
+        return self._difference_vector(
+            p1, p2, coords=coords, periodic=[False] * self.dim, axes_bounds=None
+        )
+
+    def difference_vector_real(self, p1: np.ndarray, p2: np.ndarray) -> np.ndarray:
+        # deprecated on 2024-01-09
+        warnings.warn(
+            "`difference_vector_real` has been renamed to `difference_vector`",
+            DeprecationWarning,
+        )
+        return self.difference_vector(p1, p2)
+
+    def distance(
+        self, p1: np.ndarray, p2: np.ndarray, *, coords: CoordsType = "grid"
+    ) -> float:
         """Calculate the distance between two points given in real coordinates
 
         This takes periodic boundary conditions into account if necessary.
 
-        Warning:
-            This function calculates the Euclidean distance and might thus give wrong
-            results for coordinates given in curvilinear coordinate systems.
-
         Args:
             p1 (:class:`~numpy.ndarray`):
                 First position
             p2 (:class:`~numpy.ndarray`):
                 Second position
+            coords (str):
+                The coordinate system in which the points are specified. Valid values are
+                `cartesian`, `cell`, and `grid`; see :meth:`~pde.grids.base.GridBase.transform`.
 
         Returns:
             float: Distance between the two positions
         """
-        diff = self.difference_vector_real(p1, p2)
+        diff = self.difference_vector(p1, p2, coords=coords)
         return np.linalg.norm(diff, axis=-1)  # type: ignore
 
+    def distance_real(self, p1: np.ndarray, p2: np.ndarray) -> float:
+        # deprecated on 2024-01-09
+        warnings.warn(
+            "`distance_real` has been renamed to `distance`",
+            DeprecationWarning,
+        )
+        return self.distance(p1, p2)
+
     def _get_basis(
         self, points: np.ndarray, *, coords: CoordsType = "grid"
     ) -> np.ndarray:
@@ -569,7 +625,7 @@ def _get_basis(
         Args:
             points (:class:`~numpy.ndarray`):
                 Coordinates of the point(s) where the basis determined
-            coords (:class:`~numpy.ndarray`):
+            coords (str):
                 Coordinate system in which points are specified. Valid values are
                 `cartesian`, `grid`, and `cell`; see :meth:`~pde.grids.base.GridBase.transform`.
 
@@ -666,7 +722,7 @@ def _vector_to_cartesian(
                 The coordinates of the point(s) where the vectors are specified.
             components (:class:`~numpy.ndarray`):
                 The components of the vectors at the given points
-            coords (:class:`~numpy.ndarray`):
+            coords (str):
                 The coordinate system in which the point is specified. Valid values are
                 `cartesian`, `cell`, and `grid`; see :meth:`~pde.grids.base.GridBase.transform`.
 
@@ -753,38 +809,6 @@ def normalize_point(
 
         return point
 
-    def _grid_to_cell(
-        self, grid_coords: np.ndarray, *, truncate: bool = True, normalize: bool = False
-    ) -> np.ndarray:
-        """convert grid coordinates to cell coordinates
-
-        Args:
-            grid_coords (:class:`~numpy.ndarray`):
-                The grid coordinates to convert
-            truncate (bool):
-                Flag indicating whether the resulting cell coordinates are integers
-                marking what cell the point belongs to or whether fractional coordinates
-                are returned. The default is to return integers.
-            normalize (bool):
-                Flag indicating whether the points should be normalized, e.g., whether
-                periodic boundary conditions should be imposed.
-
-        Returns:
-            :class:`~numpy.ndarray`: The cell coordinates
-        """
-        if normalize:
-            coords = self.normalize_point(grid_coords)
-        else:
-            coords = np.atleast_1d(grid_coords)
-
-        # convert from grid coordinates to cells indices
-        c_min = np.array(self.axes_bounds)[:, 0]
-        cells = (coords - c_min) / self.discretization
-        if truncate:
-            return cells.astype(int)  # type: ignore
-        else:
-            return cells  # type: ignore
-
     def transform(
         self,
         coordinates: np.ndarray,
@@ -850,7 +874,8 @@ def transform(
             if target == "grid":
                 return grid_coords
             if target == "cell":
-                return self._grid_to_cell(grid_coords, truncate=False, normalize=False)
+                c_min = np.array(self.axes_bounds)[:, 0]
+                return (grid_coords - c_min) / self.discretization  # type: ignore
 
         elif source == "cell":
             # Cell coordinates given
@@ -881,7 +906,8 @@ def transform(
             if target == "cartesian":
                 return self.point_to_cartesian(grid_coords, full=full)
             elif target == "cell":
-                return self._grid_to_cell(grid_coords, truncate=False, normalize=False)
+                c_min = np.array(self.axes_bounds)[:, 0]
+                return (grid_coords - c_min) / self.discretization  # type: ignore
             elif target == "grid":
                 return grid_coords
 
@@ -889,7 +915,6 @@ def transform(
             raise ValueError(f"Unknown source coordinates `{source}`")
         raise ValueError(f"Unknown target coordinates `{target}`")
 
-    @abstractmethod
     def polar_coordinates_real(
         self, origin: np.ndarray, *, ret_angle: bool = False
     ) -> np.ndarray | tuple[np.ndarray, ...]:
@@ -903,6 +928,7 @@ def polar_coordinates_real(
                 `False` only the distance to the origin is returned for each support
                 point of the grid. If `True`, the distance and angles are returned.
         """
+        raise NotImplementedError
 
     def contains_point(
         self,
diff --git a/pde/grids/cartesian.py b/pde/grids/cartesian.py
index 4e94c2b0..6566dfb0 100644
--- a/pde/grids/cartesian.py
+++ b/pde/grids/cartesian.py
@@ -7,6 +7,7 @@
 from __future__ import annotations
 
 import itertools
+import warnings
 from typing import TYPE_CHECKING, Any, Generator, Sequence
 
 import numpy as np
@@ -283,6 +284,13 @@ def get_random_point(
         else:
             raise ValueError(f"Unknown coordinate system `{coords}`")
 
+    def difference_vector(
+        self, p1: np.ndarray, p2: np.ndarray, *, coords: CoordsType = "grid"
+    ) -> np.ndarray:
+        return self._difference_vector(
+            p1, p2, coords=coords, periodic=self.periodic, axes_bounds=self.axes_bounds
+        )
+
     def get_line_data(self, data: np.ndarray, extract: str = "auto") -> dict[str, Any]:
         """return a line cut through the given data
 
@@ -461,12 +469,16 @@ def polar_coordinates_real(
             :class:`~numpy.ndarray` or tuple of :class:`~numpy.ndarray`:
                 Coordinates values in polar coordinates
         """
+        # deprecated on 2024-01-09
+        warnings.warn(
+            "`polar_coordinates_real` will be removed soon", DeprecationWarning
+        )
         origin = np.array(origin, dtype=np.double, ndmin=1)
         if len(origin) != self.dim:
             raise DimensionError("Dimensions are not compatible")
 
         # calculate the difference vector between all cells and the origin
-        diff = self.difference_vector_real(origin, self.cell_coords)
+        diff = self.difference_vector(origin, self.cell_coords)
         dist: np.ndarray = np.linalg.norm(diff, axis=-1)  # get distance
 
         # determine distance and optionally angles for these vectors
diff --git a/pde/grids/cylindrical.py b/pde/grids/cylindrical.py
index fc55da22..c8b542ab 100644
--- a/pde/grids/cylindrical.py
+++ b/pde/grids/cylindrical.py
@@ -6,6 +6,7 @@
 
 from __future__ import annotations
 
+import warnings
 from typing import TYPE_CHECKING, Any, Generator, Literal, Sequence
 
 import numpy as np
@@ -268,6 +269,13 @@ def get_random_point(
         else:
             raise ValueError(f"Unknown coordinate system `{coords}`")
 
+    def difference_vector(
+        self, p1: np.ndarray, p2: np.ndarray, *, coords: CoordsType = "grid"
+    ) -> np.ndarray:
+        return self._difference_vector(
+            p1, p2, coords=coords, periodic=self.periodic, axes_bounds=self.axes_bounds
+        )
+
     def get_line_data(self, data: np.ndarray, extract: str = "auto") -> dict[str, Any]:
         """return a line cut for the cylindrical grid
 
@@ -459,6 +467,10 @@ def polar_coordinates_real(
                 each support point of the grid. If `True`, the distance and angles are
                 returned.
         """
+        # deprecated on 2024-01-09
+        warnings.warn(
+            "`polar_coordinates_real` will be removed soon", DeprecationWarning
+        )
         origin = np.array(origin, dtype=np.double, ndmin=1)
         if len(origin) != self.dim:
             raise DimensionError("Dimensions are not compatible")
@@ -467,7 +479,7 @@ def polar_coordinates_real(
             raise RuntimeError("Origin must lie on symmetry axis for cylindrical grid")
 
         # calculate the difference vector between all cells and the origin
-        diff = self.difference_vector_real(np.array([0, origin[2]]), self.cell_coords)
+        diff = self.difference_vector(np.array([0, origin[2]]), self.cell_coords)
         dist: np.ndarray = np.linalg.norm(diff, axis=-1)  # get distance
 
         if ret_angle:
diff --git a/pde/grids/spherical.py b/pde/grids/spherical.py
index aa9863e9..e65e0751 100644
--- a/pde/grids/spherical.py
+++ b/pde/grids/spherical.py
@@ -10,6 +10,7 @@
 
 from __future__ import annotations
 
+import warnings
 from abc import ABCMeta
 from typing import TYPE_CHECKING, Any, Literal, TypeVar
 
@@ -371,6 +372,10 @@ def polar_coordinates_real(
                 point of the grid. If `True`, the distance and angles are returned. Note
                 that in the case of spherical grids, this angle is zero by convention.
         """
+        # deprecated on 2024-01-09
+        warnings.warn(
+            "`polar_coordinates_real` will be removed soon", DeprecationWarning
+        )
         # check the consistency of the origin argument, which can be set for other grids
         if origin is not None:
             origin = np.array(origin, dtype=np.double, ndmin=1)
diff --git a/tests/grids/test_cartesian_grids.py b/tests/grids/test_cartesian_grids.py
index cf5f0b51..70a30fb8 100644
--- a/tests/grids/test_cartesian_grids.py
+++ b/tests/grids/test_cartesian_grids.py
@@ -100,7 +100,7 @@ def norm_numba_wrap(x):
     p = rng.random(1) * grid.shape
     d, a = grid.polar_coordinates_real(p, ret_angle=True)
     c2 = grid.from_polar_coordinates(d, a, p)
-    assert np.allclose(grid.distance_real(c1, c2), 0)
+    assert np.allclose(grid.distance(c1, c2), 0)
 
     # test boundary points
     np.testing.assert_equal(grid._boundary_coordinates(0, False), np.array([0]))
@@ -138,7 +138,7 @@ def test_unit_grid_2d(rng):
     p = rng.random(2) * grid.shape
     d, a = grid.polar_coordinates_real(p, ret_angle=True)
     c2 = grid.from_polar_coordinates(d, a, p)
-    assert np.allclose(grid.distance_real(c1, c2), 0)
+    assert np.allclose(grid.distance(c1, c2), 0)
 
     # test boundary points
     np.testing.assert_equal(
@@ -221,7 +221,7 @@ def test_rect_grid_1d(rng):
         p = rng.random(1) * grid.shape
         d, a = grid.polar_coordinates_real(p, ret_angle=True)
         c2 = grid.from_polar_coordinates(d, a, p)
-        assert np.allclose(grid.distance_real(c1, c2), 0)
+        assert np.allclose(grid.distance(c1, c2), 0)
 
 
 def test_rect_grid_2d(rng):
@@ -252,7 +252,7 @@ def test_rect_grid_2d(rng):
     d, a = grid.polar_coordinates_real(p, ret_angle=True)
     c2 = grid.from_polar_coordinates(d, a, p)
 
-    assert np.allclose(grid.distance_real(c1, c2), 0)
+    assert np.allclose(grid.distance(c1, c2), 0)
 
 
 def test_rect_grid_3d(rng):
@@ -297,7 +297,7 @@ def test_unit_rect_grid(periodic, rng):
 
     for _ in range(10):
         p1, p2 = rng.normal(scale=10, size=(2, dim))
-        assert g1.distance_real(p1, p2) == pytest.approx(g2.distance_real(p1, p2))
+        assert g1.distance(p1, p2) == pytest.approx(g2.distance(p1, p2))
 
     p0 = rng.normal(scale=10, size=dim)
     np.testing.assert_allclose(