diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 4d0f9c8d954c8376dc40aaf128fcde879fbc02a1..d91600c751c3d81b42c9f2c4cea29185c1a9c34f 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -170,26 +170,28 @@ def separate_into_bulk_and_interface(free_energy):
 
 
 def analytic_interface_profile(x, interface_width=interface_width_symbol):
-    """Analytic expression for a 1D interface normal to x with given interface width
+    """Analytic expression for a 1D interface normal to x with given interface width.
 
     The following doctest shows that the returned analytical solution is indeed a solution of the ODE that we
     get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
     parameter, i.e. at a transition between two phases.
-    >>> numPhases = 4
-    >>> x, phi = sp.Symbol("x"), symbolic_order_parameters(numPhases-1)
-    >>> F = free_energy_functional_n_phases(order_parameters=phi)
-    >>> mu = chemical_potentials_from_free_energy(F)
-    >>> mu0 = mu[0].subs({p: 0 for p in phi[1:]})  # mu[0] as function of one order parameter only
-    >>> solution = analytic_interface_profile(x)
-    >>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
-    >>> sp.expand(mu0.subs(solutionSubstitution))  # inserting solution should solve the mu_0=0 equation
-    0
+
+    Examples:
+        >>> numPhases = 4
+        >>> x, phi = sp.Symbol("x"), symbolic_order_parameters(numPhases-1)
+        >>> F = free_energy_functional_n_phases(order_parameters=phi)
+        >>> mu = chemical_potentials_from_free_energy(F)
+        >>> mu0 = mu[0].subs({p: 0 for p in phi[1:]})  # mu[0] as function of one order parameter only
+        >>> solution = analytic_interface_profile(x)
+        >>> solutionSubstitution = {phi[0]: solution, Diff(Diff(phi[0])): sp.diff(solution, x, x) }
+        >>> sp.expand(mu0.subs(solutionSubstitution))  # inserting solution should solve the mu_0=0 equation
+        0
     """
     return (1 + sp.tanh(x / (2 * interface_width))) / 2
 
 
 def chemical_potentials_from_free_energy(free_energy, order_parameters=None):
-    """Computes chemical potentials as functional derivative of free energy"""
+    """Computes chemical potentials as functional derivative of free energy."""
     symbols = free_energy.atoms(sp.Symbol)
     if order_parameters is None:
         order_parameters = [s for s in symbols if s.name.startswith(orderParameterSymbolName)]
diff --git a/phasefield/experiments1D.py b/phasefield/experiments1D.py
index b2360cc270f81d108050aa1db1b35734f88758ab..110bb8cadb0bed3b0acbf4dfd76d8d723456c1c9 100644
--- a/phasefield/experiments1D.py
+++ b/phasefield/experiments1D.py
@@ -13,26 +13,26 @@ def plot_status(phase_field_step, from_x=None, to_x=None):
 
     dh = phase_field_step.data_handling
 
-    num_phases = phase_field_step.numOrderParameters
+    num_phases = phase_field_step.num_order_parameters
 
     plt.subplot(1, 3, 1)
     plt.title('φ')
-    phi_name = phase_field_step.phiFieldName
+    phi_name = phase_field_step.phi_field_name
     for i in range(num_phases):
         plt.plot(dh.gather_array(phi_name, make_slice[from_x:to_x, 0, i]), marker='x', label='φ_%d' % (i,))
     plt.legend()
 
     plt.subplot(1, 3, 2)
     plt.title("μ")
-    mu_name = phase_field_step.muFieldName
+    mu_name = phase_field_step.mu_field_name
     for i in range(num_phases):
         plt.plot(dh.gather_array(mu_name, make_slice[from_x:to_x, 0, i]), marker='x', label='μ_%d' % (i,))
     plt.legend()
 
     plt.subplot(1, 3, 3)
     plt.title("Force and Velocity")
-    plt.plot(dh.gather_array(phase_field_step.forceFieldName, make_slice[from_x:to_x, 0, 0]), label='F', marker='x')
-    plt.plot(dh.gather_array(phase_field_step.velFieldName, make_slice[from_x:to_x, 0, 0]), label='u', marker='v')
+    plt.plot(dh.gather_array(phase_field_step.force_field_name, make_slice[from_x:to_x, 0, 0]), label='F', marker='x')
+    plt.plot(dh.gather_array(phase_field_step.vel_field_name, make_slice[from_x:to_x, 0, 0]), label='u', marker='v')
     plt.legend()
 
 
@@ -59,14 +59,14 @@ def init_sharp_interface(pf_step, phase_idx, x1=None, x2=None, inverse=False):
     if x2 is None:
         x2 = 3 * x1
 
-    if phase_idx >= pf_step.numOrderParameters:
+    if phase_idx >= pf_step.num_order_parameters:
         return
 
     for b in pf_step.data_handling.iterate():
         x = b.cell_index_arrays[0]
         mid = np.logical_and(x1 < x, x < x2)
 
-        phi = b[pf_step.phiFieldName]
+        phi = b[pf_step.phi_field_name]
         val1, val2 = (1, 0) if inverse else (0, 1)
 
         phi[..., phase_idx].fill(val1)
@@ -94,7 +94,7 @@ def tanh_test(pf_step, phase0, phase1, expected_interface_width=1, time_steps=10
 
     domain_size = pf_step.data_handling.shape
     pf_step.reset()
-    pf_step.data_handling.fill(pf_step.phiFieldName, 0)
+    pf_step.data_handling.fill(pf_step.phi_field_name, 0)
     init_sharp_interface(pf_step, phase_idx=phase0, inverse=False)
     init_sharp_interface(pf_step, phase_idx=phase1, inverse=True)
     pf_step.set_pdf_fields_from_macroscopic_values()
@@ -140,14 +140,14 @@ def galilean_invariance_test(pf_step, velocity=0.05, rounds=3, phase0=0, phase1=
     print("Velocity:", velocity, " Time steps for round:", round_time_steps)
 
     pf_step.reset()
-    pf_step.data_handling.fill(pf_step.phiFieldName, 0)
+    pf_step.data_handling.fill(pf_step.phi_field_name, 0)
     init_sharp_interface(pf_step, phase_idx=phase0, inverse=False)
     init_sharp_interface(pf_step, phase_idx=phase1, inverse=True)
     pf_step.set_pdf_fields_from_macroscopic_values()
 
     print("Running", init_time_steps, "initial time steps")
     pf_step.run(init_time_steps)
-    pf_step.data_handling.fill(pf_step.velFieldName, velocity, value_idx=0)
+    pf_step.data_handling.fill(pf_step.vel_field_name, velocity, value_idx=0)
     pf_step.set_pdf_fields_from_macroscopic_values()
 
     step_location = domain_size[0] // 4
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index a787062a4c999db3eaca708d8a877b00d634ede0..b32c70e4df4d4d48c76710d06e0c223a8e427010 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -30,11 +30,11 @@ class PhaseFieldStep:
         if data_handling is None:
             data_handling = SerialDataHandling(domain_size, periodicity=True)
 
-        self.freeEnergy = free_energy
-        self.concentrationToOrderParameter = concentration_to_order_parameters
-        self.orderParametersToConcentrations = order_parameters_to_concentrations
+        self.free_energy = free_energy
+        self.concentration_to_order_parameter = concentration_to_order_parameters
+        self.order_parameters_to_concentrations = order_parameters_to_concentrations
 
-        self.chemicalPotentials = chemical_potentials_from_free_energy(free_energy, order_parameters)
+        self.chemical_potentials = chemical_potentials_from_free_energy(free_energy, order_parameters)
 
         # ------------------ Adding arrays ---------------------
         gpu = target == 'gpu'
@@ -61,7 +61,7 @@ class PhaseFieldStep:
 
         # ------------------ Creating kernels ------------------
         phi = tuple(self.phi_field(i) for i in range(len(order_parameters)))
-        F = self.freeEnergy.subs({old: new for old, new in zip(order_parameters, phi)})
+        F = self.free_energy.subs({old: new for old, new in zip(order_parameters, phi)})
 
         if homogeneous_neumann_boundaries:
             def apply_neumann_boundaries(eqs):
@@ -77,7 +77,7 @@ class PhaseFieldStep:
         # μ and pressure tensor update
         self.phiSync = data_handling.synchronization_function([self.phi_field_name], target=target)
         self.muEqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx)
-        self.pressureTensorEqs = pressure_tensor_kernel(self.freeEnergy, order_parameters,
+        self.pressureTensorEqs = pressure_tensor_kernel(self.free_energy, order_parameters,
                                                         self.phi_field, self.pressure_tensor_field, dx)
         mu_and_pressure_tensor_eqs = self.muEqs + self.pressureTensorEqs
         mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs)
@@ -225,8 +225,8 @@ class PhaseFieldStep:
         return self.hydroLbmStep.boundary_handling
 
     def set_concentration(self, slice_obj, concentration):
-        if self.concentrationToOrderParameter is not None:
-            phi = self.concentrationToOrderParameter(concentration)
+        if self.concentration_to_order_parameter is not None:
+            phi = self.concentration_to_order_parameter(concentration)
         else:
             phi = np.array(concentration)
 
@@ -255,7 +255,7 @@ class PhaseFieldStep:
 
     def concentration_slice(self, slice_obj=None):
         phi = self.phi_slice(slice_obj)
-        return phi if self.orderParametersToConcentrations is None else self.orderParametersToConcentrations(phi)
+        return phi if self.order_parameters_to_concentrations is None else self.order_parameters_to_concentrations(phi)
 
     def mu_slice(self, slice_obj=None):
         return self._get_slice(self.mu_field_name, slice_obj)