diff --git a/PySDM/formulae.py b/PySDM/formulae.py
index 4350ab40e..b5fdc247e 100644
--- a/PySDM/formulae.py
+++ b/PySDM/formulae.py
@@ -83,7 +83,9 @@ def __init__(  # pylint: disable=too-many-locals
         self.air_dynamic_viscosity = air_dynamic_viscosity
 
         self._components = tuple(
-            i for i in dir(self) if not i.startswith("__") and i != "flatten"
+            i
+            for i in dir(self)
+            if not i.startswith("__") and i not in ("flatten", "get_constant")
         )
 
         constants_defaults = {
@@ -161,6 +163,11 @@ def flatten(self):
             functions[attr] = getattr(self, attr)
         return namedtuple("FlattenedFormulae", functions.keys())(**functions)
 
+    def get_constant(self, key: str):
+        """getter-like method for cases where using the `constants` named tuple is not possible
+        (e.g., if calling from a language which does not support named tuples)"""
+        return getattr(self.constants, key)
+
 
 def _formula(func, constants, dimensional_analysis, **kw):
     parameters_keys = tuple(inspect.signature(func).parameters.keys())
diff --git a/README.md b/README.md
index b7951c188..c97d3989a 100644
--- a/README.md
+++ b/README.md
@@ -321,14 +321,15 @@ In the listing below, its usage is interleaved with plotting logic
 <summary>Julia (click to expand)</summary>
 
 ```Julia
-rho_w = pyimport("PySDM.physics.constants_defaults").rho_w
 using Plots; plotlyjs()
 
 for step = 0:1200:3600
     particulator.run(step - particulator.n_steps)
     plot!(
         radius_bins_edges[1:end-1] / si.um,
-        particulator.products["dv/dlnr"].get()[:] * rho_w / si.g,
+        particulator.formulae.particle_shape_and_density.volume_to_mass(
+            particulator.products["dv/dlnr"].get()[:]
+        )/ si.g,
         linetype=:steppost,
         xaxis=:log,
         xlabel="particle radius [µm]",
@@ -343,12 +344,12 @@ savefig("plot.svg")
 <summary>Matlab (click to expand)</summary>
 
 ```Matlab
-rho_w = py.importlib.import_module('PySDM.physics.constants_defaults').rho_w;
-
 for step = 0:1200:3600
     particulator.run(int32(step - particulator.n_steps));
     x = radius_bins_edges / si.um;
-    y = particulator.products{"dv/dlnr"}.get() * rho_w / si.g;
+    y = particulator.formulae.particle_shape_and_density.volume_to_mass( ...
+        particulator.products{"dv/dlnr"}.get() ...
+    ) / si.g;
     stairs(...
         x(1:end-1), ... 
         double(py.array.array('d',py.numpy.nditer(y))), ...
@@ -367,14 +368,17 @@ legend()
 <summary>Python (click to expand)</summary>
 
 ```Python
-from PySDM.physics.constants_defaults import rho_w
 from matplotlib import pyplot
 
 for step in [0, 1200, 2400, 3600]:
     particulator.run(step - particulator.n_steps)
-    pyplot.step(x=radius_bins_edges[:-1] / si.um,
-                y=particulator.products['dv/dlnr'].get()[0] * rho_w / si.g,
-                where='post', label=f"t = {step}s")
+    pyplot.step(
+        x=radius_bins_edges[:-1] / si.um,
+        y=particulator.formulae.particle_shape_and_density.volume_to_mass(
+            particulator.products['dv/dlnr'].get()[0]
+        ) / si.g,
+        where='post', label=f"t = {step}s"
+    )
 
 pyplot.xscale('log')
 pyplot.xlabel('particle radius [µm]')
diff --git a/tests/unit_tests/test_formulae.py b/tests/unit_tests/test_formulae.py
index f5e80c6bc..b146cf787 100644
--- a/tests/unit_tests/test_formulae.py
+++ b/tests/unit_tests/test_formulae.py
@@ -78,3 +78,14 @@ def test_flatten():
         # assert
         temp = 300 * si.K
         assert sut.latent_heat__lv(temp) == f.latent_heat.lv(temp)
+
+    @staticmethod
+    def test_get_constant():
+        # arrange
+        rho_w = 666 * si.kg / si.m**3
+
+        # act
+        sut = formulae.Formulae(constants={"rho_w": rho_w})
+
+        # assert
+        assert sut.get_constant("rho_w") == rho_w