diff --git a/src/rt1_model/_calc.py b/src/rt1_model/_calc.py
index f4b840c..a2cc6c9 100644
--- a/src/rt1_model/_calc.py
+++ b/src/rt1_model/_calc.py
@@ -564,7 +564,11 @@ def calc(self, **params):
                 {
                     *self._all_param_symbs,
                     *map(
-                        str, [*self.V._func.free_symbols, *self.SRF._func.free_symbols]
+                        str,
+                        [
+                            *self.V.phase_function.free_symbols,
+                            *self.SRF.phase_function.free_symbols,
+                        ],
                     ),
                 }
                 - {"phi_0", "phi_ex", "theta_0", "theta_ex"}
@@ -1415,7 +1419,7 @@ def _d_volume_dummy_lambda(self, key):
 
         return _lambdify(
             args,
-            sp.diff(self.V._func, sp.Symbol(key)),
+            sp.diff(self.V.phase_function, sp.Symbol(key)),
         )
 
     def _d_surface_ddummy(self, key):
@@ -1486,7 +1490,7 @@ def _d_surface_dummy_lambda(self, key):
 
         return _lambdify(
             args,
-            sp.diff(self.SRF._func, sp.Symbol(key)),
+            sp.diff(self.SRF.phase_function, sp.Symbol(key)),
         )
 
     def jacobian(self, param_list=["omega", "tau", "NormBRDF"], format="list"):
diff --git a/src/rt1_model/_scatter.py b/src/rt1_model/_scatter.py
index 0793265..3c10384 100644
--- a/src/rt1_model/_scatter.py
+++ b/src/rt1_model/_scatter.py
@@ -90,16 +90,16 @@ def _lambda_func(self, *args):
 
         # replace arguments and evaluate expression
         args = (theta_0, theta_ex, phi_0, phi_ex) + tuple(args)
-        pfunc = _lambdify(args, self._func)
+        pfunc = _lambdify(args, self.phase_function)
 
         # TODO check requirement for this!
-        # if _func is a constant, lambdify will create a function that returns a scalar
+        # if phase_function is a constant, lambdify will create a function that returns a scalar
         # which is not suitable for further processing. in that case, vectorize the
         # obtained function
 
         # TODO maybe find a better check for this
-        # if self._func.is_constant():   # this is too slow!
-        # if len(self._func.free_symbols) == 0:
+        # if self.phase_function.is_constant():   # this is too slow!
+        # if len(self.phase_function.free_symbols) == 0:
         #     pfunc = np.vectorize(pfunc)
 
         return pfunc
@@ -246,13 +246,13 @@ def __init__(self, choices=None, **kwargs):
         super().__init__(**kwargs)
 
     @property
-    def _func(self):
-        """Phase function as sympy object for later evaluation."""
-        _func = 0
+    def phase_function(self):
+        """Phase function as sympy expression for later evaluation."""
+        func = 0
         for c, o in zip(self._weights, self._objs):
-            _func += c * o._func
+            func += c * o.phase_function
 
-        return _func
+        return func
 
     @property
     def ncoefs(self):
diff --git a/src/rt1_model/plot.py b/src/rt1_model/plot.py
index 7dabaef..fec9631 100644
--- a/src/rt1_model/plot.py
+++ b/src/rt1_model/plot.py
@@ -113,7 +113,7 @@ def polarplot(
         param_dict = [param_dict]
 
     # Check if all required parameters have been provided in the param_dict
-    required_symbs = set(map(str, V_SRF._func.free_symbols)) - {
+    required_symbs = set(map(str, V_SRF.phase_function.free_symbols)) - {
         "phi_0",
         "phi_ex",
         "theta_0",
@@ -397,7 +397,9 @@ def _check_params(R, param_dict, additional_params=[]):
     # check if all required parameters for the analyzers have been defined
     symbs = {
         *R._all_param_symbs,
-        *map(str, [*R.V._func.free_symbols, *R.SRF._func.free_symbols]),
+        *map(
+            str, [*R.V.phase_function.free_symbols, *R.SRF.phase_function.free_symbols]
+        ),
     } - {"phi_0", "phi_ex", "theta_0", "theta_ex"}
     for p in additional_params:
         symbs.add(p)
diff --git a/src/rt1_model/surface.py b/src/rt1_model/surface.py
index 5865533..4545dc5 100644
--- a/src/rt1_model/surface.py
+++ b/src/rt1_model/surface.py
@@ -48,7 +48,7 @@ def legendre_coefficients(self):
         ...
 
     @abstractmethod
-    def _func(self):
+    def phase_function(self):
         """Phase function as sympy object."""
         ...
 
@@ -93,7 +93,7 @@ def legendre_coefficients(self):
         return (1.0 / sp.pi) * sp.KroneckerDelta(0, n)
 
     @property
-    def _func(self):
+    def phase_function(self):
         """Phase function as sympy object."""
         return 1.0 / sp.pi
 
@@ -141,7 +141,7 @@ def legendre_coefficients(self):
         )
 
     @property
-    def _func(self):
+    def phase_function(self):
         """Phase function as sympy object."""
         return (
             1.0
@@ -178,7 +178,7 @@ def __init__(self, t=None, **kwargs):
         self.t = _parse_sympy_param(t)
 
     @property
-    def _func(self):
+    def phase_function(self):
         """Phase function as sympy object."""
 
         return (
@@ -221,7 +221,7 @@ def __init__(self, t=None, **kwargs):
         self.t = _parse_sympy_param(t)
 
     @property
-    def _func(self):
+    def phase_function(self):
         """Define Phase function as sympy object."""
 
         nadir_hemreflect = 4 * (
diff --git a/src/rt1_model/volume.py b/src/rt1_model/volume.py
index 05a6380..14af92b 100644
--- a/src/rt1_model/volume.py
+++ b/src/rt1_model/volume.py
@@ -47,8 +47,8 @@ def legendre_coefficients(self):
         ...
 
     @abstractmethod
-    def _func(self):
-        """Phase function as sympy object."""
+    def phase_function(self):
+        """Phase function as sympy expression."""
         ...
 
 
@@ -92,8 +92,8 @@ def legendre_coefficients(self):
         return (1.0 / (4.0 * sp.pi)) * sp.KroneckerDelta(0, n)
 
     @property
-    def _func(self):
-        """Phase function as sympy object."""
+    def phase_function(self):
+        """Phase function as sympy expression."""
         return 1.0 / (4.0 * sp.pi)
 
 
@@ -118,8 +118,8 @@ def ncoefs(self):
         return 3
 
     @property
-    def _func(self):
-        """Phase function as sympy object."""
+    def phase_function(self):
+        """Phase function as sympy expression."""
         return 3.0 / (16.0 * sp.pi) * (1.0 + self.scat_angle_symbolic**2.0)
 
     @property
@@ -160,8 +160,8 @@ def __init__(self, t=None, **kwargs):
         self.t = _parse_sympy_param(t)
 
     @property
-    def _func(self):
-        """Phase function as sympy object."""
+    def phase_function(self):
+        """Phase function as sympy expression."""
         func = (1.0 - self.t**2.0) / (
             (4.0 * sp.pi)
             * (1.0 + self.t**2.0 - 2.0 * self.t * self.scat_angle_symbolic) ** 1.5
@@ -201,8 +201,8 @@ def __init__(self, t=None, **kwargs):
         self.t = _parse_sympy_param(t)
 
     @property
-    def _func(self):
-        """Phase function as sympy object."""
+    def phase_function(self):
+        """Phase function as sympy expression."""
         return (
             3.0
             / (8.0 * sp.pi)
diff --git a/tests/test_volume_surface.py b/tests/test_volume_surface.py
index 76643d0..383385b 100644
--- a/tests/test_volume_surface.py
+++ b/tests/test_volume_surface.py
@@ -53,7 +53,7 @@ def test_surface_init(backend):
 
         SRF.calc(0.1, 0.2, 0.3, 0.4)
         SRF.legendre_expansion(0.1, 0.2, 0.3, 0.4)
-        SRF._func
+        SRF.phase_function
 
 
 @pytest.mark.parametrize("backend", ["sympy", "symengine"])
@@ -85,7 +85,7 @@ def test_volume_init(backend):
         # evaluate function numerical
         V.calc(0.1, 0.2, 0.3, 0.4)
         V.legendre_expansion(0.1, 0.2, 0.3, 0.4)
-        V._func
+        V.phase_function
 
 
 @pytest.mark.parametrize("backend", ["sympy", "symengine"])
@@ -110,7 +110,7 @@ def test_linear_combinations_SRF(backend):
 
     SRF.calc(0.1, 0.2, 0.3, 0.4)
     SRF.legendre_expansion(0.1, 0.2, 0.3, 0.4)
-    SRF._func
+    SRF.phase_function
 
 
 @pytest.mark.parametrize("backend", ["sympy", "symengine"])
@@ -135,4 +135,4 @@ def test_linear_combinations_V(backend):
 
     V.calc(0.1, 0.2, 0.3, 0.4)
     V.legendre_expansion(0.1, 0.2, 0.3, 0.4)
-    V._func
+    V.phase_function