diff --git a/phoebe/backend/backends.py b/phoebe/backend/backends.py
index 1b9d10873..352e49d45 100644
--- a/phoebe/backend/backends.py
+++ b/phoebe/backend/backends.py
@@ -14,7 +14,7 @@
 from phoebe.parameters import StringParameter, DictParameter, ArrayParameter, ParameterSet
 from phoebe.parameters.parameters import _extract_index_from_string
 from phoebe import dynamics
-from phoebe.backend import universe, horizon_analytic
+from phoebe.backend import universe, horizon_analytic, contacts_energy_transfer
 from phoebe.atmospheres import passbands
 from phoebe.distortions  import roche
 from phoebe.frontend import io
@@ -976,6 +976,13 @@ def _compute_intrinsic_system_at_t0(self, b, compute,
         # kinds = [b.get_dataset(dataset=ds).kind for ds in datasets]
 
         system.update_positions(t0, x0, y0, z0, vx0, vy0, vz0, etheta0, elongan0, eincl0, ignore_effects=True)
+
+        #ENERGY TRANSFER FOR CONTACTS
+        if 'envelope' in b.filter(context='component'):  # only makes sense when an envelope is present
+            mixing_enabled = b.get_value(qualifier='mixing_enabled', context='component', **_skip_filter_checks)
+            if mixing_enabled:
+                self._do_mixing(b, system)
+
         system.populate_observables(t0, ['lc' for dataset in datasets], datasets, ignore_effects=True)
 
 
@@ -985,6 +992,33 @@ def _compute_intrinsic_system_at_t0(self, b, compute,
 
         return system
 
+    @staticmethod
+    def _do_mixing(b, system):
+        mixing_method = b.get_value(qualifier='mixing_method', context='component', **_skip_filter_checks)
+        mixing_power = b.get_value(qualifier='mixing_power', context='component', **_skip_filter_checks)
+        secondary_teff = b.get_value(qualifier='teff', context='component', component='secondary',
+                                     **_skip_filter_checks)
+        primary_teff = b.get_value(qualifier='teff', context='component', component='primary', **_skip_filter_checks)
+        teff_ratio = secondary_teff / primary_teff
+
+        primary_mesh, secondary_mesh = system.bodies[0].meshes.values()
+        coords1 = primary_mesh.roche_coords_for_computations
+        teffs1 = primary_mesh.teffs
+        coords2 = secondary_mesh.roche_coords_for_computations
+        teffs2 = secondary_mesh.teffs
+
+        new_teffs1, new_teffs2 = contacts_energy_transfer.mix_teffs(np.array(coords1), np.array(teffs1),
+                                                                    np.array(coords2), np.array(teffs2),
+                                                                    mixing_method=mixing_method,
+                                                                    mixing_power=mixing_power, teff_ratio=teff_ratio)
+        # w=smoothing_factor, cutoff=0.)
+        primary_mesh.update_columns(teffs=new_teffs1)
+        secondary_mesh.update_columns(teffs=new_teffs2)
+
+        # print("new_teffs1.median", np.median(new_teffs1))
+        # print("new_teffs2.median", np.median(new_teffs2))
+        system.bodies[0]._halves[0].smoothed_teffs = new_teffs1
+        system.bodies[0]._halves[1].smoothed_teffs = new_teffs2
 
     def _worker_setup(self, b, compute, times, infolists, **kwargs):
         logger.debug("rank:{}/{} PhoebeBackend._worker_setup: extracting parameters".format(mpi.myrank, mpi.nprocs))
@@ -1072,6 +1106,12 @@ def _run_single_time(self, b, i, time, infolist, **kwargs):
             logger.debug("rank:{}/{} PhoebeBackend._run_single_time: calling system.update_positions at time={}".format(mpi.myrank, mpi.nprocs, time))
             system.update_positions(time, xi, yi, zi, vxi, vyi, vzi, ethetai, elongani, eincli, ds=di, Fs=Fi)
 
+            #ENERGY TRANSFER FOR CONTACTS
+            if 'envelope' in b.filter(context='component'):  # only makes sense when an envelope is present
+                mixing_enabled = b.get_value(qualifier='mixing_enabled', context='component', **_skip_filter_checks)
+                if mixing_enabled and i==0:
+                    self._do_mixing(b, system)
+
             # Now we need to determine which triangles are visible and handle subdivision
             # NOTE: this should come after populate_observables so that each subdivided triangle
             # will have identical local quantities.  The only downside to this is that we can't
diff --git a/phoebe/backend/contacts_energy_transfer.py b/phoebe/backend/contacts_energy_transfer.py
new file mode 100644
index 000000000..f441c3ef6
--- /dev/null
+++ b/phoebe/backend/contacts_energy_transfer.py
@@ -0,0 +1,358 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from itertools import combinations
+from scipy.optimize import curve_fit
+import logging
+
+logger = logging.getLogger("CONTACT ENVELOPES")
+logger.addHandler(logging.NullHandler())
+
+
+def _isolate_neck(coords_all, teffs_all, cutoff=0., component=1, plot=False):
+    """
+    Selects vertices in the x = [0 + cutoff, 1 - cutoff] region of the Roche geometry.
+
+    Parameters
+    ----------
+    coords_all: 3D points
+    teffs_all: Teffs to filter on
+    cutoff: wiggle room value for the edges (default=0)
+    component: which component? 1 or 2 (default=1)
+    plot: plot the resulting selection? (default=False)
+
+    Returns
+    -------
+    Filtered coordinates, filtered Teffs, indices of coords_all that was filtered on
+    """
+    if component == 1:
+        cond = coords_all[:, 0] >= 0 + cutoff
+    elif component == 2:
+        cond = coords_all[:, 0] <= 1 - cutoff
+    else:
+        raise ValueError
+
+    if plot:
+        fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
+        axes[0].scatter(coords_all[cond][:, 0], coords_all[cond][:, 1])
+        axes[1].scatter(coords_all[cond][:, 0], teffs_all[cond])
+        axes[0].set_xlabel('x')
+        axes[0].set_ylabel('y')
+        axes[1].set_xlabel('x')
+        axes[1].set_ylabel('teff')
+        fig.tight_layout()
+        plt.show()
+
+    return coords_all[cond], teffs_all[cond], np.argwhere(cond).flatten()
+
+
+def _dist(p1, p2):
+    """
+    Euclidean distance between two points in R^n
+    Parameters
+    ----------
+    p1: point 1
+    p2: point 2
+
+    Returns
+    -------
+    Square root of dot product of p2 - p1
+    """
+    return np.sqrt(np.dot(p2 - p1, p2 - p1))
+
+
+def _isolate_sigma_fitting_regions(coords_neck, teffs_neck, direction='x', cutoff=0., component=1, plot=False):
+    """
+    Determine the region on which sigma smoothing is going to be fitted. If `direction=="y"`, this is going to be the
+    slice close to `x=0` (or `x=1` for the secondary). If `direction=="x"`, this is going to be the slice near the
+    equator (`y=0`).
+
+    Parameters
+    ----------
+    coords_neck: 3D points of the neck region
+    teffs_neck: Teffs of the neck region
+    direction: which direction? either 'x' or 'y'
+    cutoff: wiggle room value for the edges (default=0)
+    component: which component? 1 or 2 (default=1)
+    plot: plot the resulting selection? (default=False)
+
+    Returns
+    -------
+    Filtered coordinates, filtered Teffs
+    """
+    distances = [_dist(p1, p2) for p1, p2 in combinations(coords_neck, 2)]
+    min_dist = np.min(distances[distances != 0])  # smallest dist. between any two vertices (~ resolution of the mesh)
+
+    if direction == 'x':
+        cond = (coords_neck[:, 1] >= -cutoff - 0.2 * min_dist) & (coords_neck[:, 1] <= cutoff + 0.2 * min_dist)
+
+    elif direction == 'y':
+        if component == 1:
+            cond = coords_neck[:, 0] <= 0 + cutoff + 0.15 * min_dist
+        elif component == 2:
+            cond = coords_neck[:, 0] >= 1 - cutoff - 0.15 * min_dist
+        else:
+            raise ValueError('Unknown component, must be 1 or 2')
+    else:
+        raise ValueError('Unknown direction, must be either "x" or "y"')
+
+    if plot:
+        if direction == 'x':
+            plt.scatter(coords_neck[cond][:, 0], teffs_neck[cond])
+            plt.show()
+        elif direction == 'y':
+            plt.scatter(coords_neck[cond][:, 1], teffs_neck[cond])
+            plt.show()
+
+    return coords_neck[cond], teffs_neck[cond]
+
+
+def _compute_new_teff_at_neck(coords1, teffs1, coords2, teffs2, w=0.5):
+    """
+    Computes the target teff at the neck that should be attained by the smoothing.
+    Parameters
+    ----------
+    coords1: Roche coordinates of the primary
+    teffs1: Teffs of the primary
+    coords2: Roche coordinates of the secondary
+    teffs2: Teffs of the secondary
+    w: Weight of the primary contribution relative to the secondary
+
+    Returns
+    -------
+    x coordinate of the neck, target Teff at the neck
+    """
+    distances1 = [_dist(p1, p2) for p1, p2 in combinations(coords1, 2)]
+    min_dist1 = np.min(distances1[distances1 != 0])
+    distances2 = [_dist(p1, p2) for p1, p2 in combinations(coords2, 2)]
+    min_dist2 = np.min(distances2[distances2 != 0])
+
+    x_neck = np.average((coords1[:, 0].max(), coords2[:, 0].min()))
+
+    teffs_neck1 = teffs1[coords1[:, 0] >= coords1[:, 0].max() - 0.25 * min_dist1]
+    teffs_neck2 = teffs2[coords2[:, 0] <= coords2[:, 0].min() + 0.25 * min_dist2]
+
+    teff1 = np.max(teffs2)  # np.average(teffs_neck1)
+    teff2 = np.average(teffs_neck2)
+    tavg = w * teff1 + (1 - w) * teff2
+    if tavg > teffs2.max():
+        logger.warning('Tavg > Teff2, setting new temperature to 1 percent of Teff2 max. %i > %i' % (
+            int(tavg), int(teffs2.max())))
+        tavg = teffs2.max() - 0.01 * teffs2.max()
+
+    return x_neck, tavg
+
+
+def _compute_sigmax(Tavg, x, x0, offset, amplitude):
+    """
+    Rather than fitting for sigma_x, we simply invert a gaussian profile to match the target Tavg
+    (minus a potential offset)
+    """
+    return ((-1) * (x - x0) ** 2 / np.log((Tavg - offset) / amplitude)) ** 0.5
+
+
+def _fit_sigma(coords, teffs, offset=0., cutoff=0., direction='y', component=1, plot=False):
+    """
+    Fits the width of a gaussian distribution to given coords and teffs.
+
+    Parameters
+    ----------
+    coords: 3D coordinates in the Roche geometry
+    teffs: teffs at the given coordinates
+    offset: temperature offset to include (default=0)
+    cutoff: x-displacement of the center of the gaussian (default=0)
+    direction: direction to fit gaussian in (default='y')
+    component: which component? 1 or 2 (default=1)
+    plot: plot the resulting computation? (default=False)
+
+    Returns
+    -------
+    the computed width, amplitude, and the fitted model Teffs
+    """
+
+    def gaussian_1d(x, s):
+        a = 1. / s ** 2
+        g = offset + amplitude * np.exp(- (a * ((x - x0) ** 2)))
+        return g
+
+    if direction == 'y':
+        coord_ind = 1
+        x0 = 0.
+    elif direction == 'x':
+        coord_ind = 0
+        if component == 1:
+            x0 = 0 + cutoff
+        elif component == 2:
+            x0 = 1 - cutoff
+        else:
+            raise ValueError('Unknown component selected, must be "1" or "2"')
+    else:
+        raise ValueError('Unknown direction selected, must be "y" or "x"')
+
+    amplitude = teffs.max() - offset
+    sigma_0 = 0.5
+    result = curve_fit(gaussian_1d, xdata=coords[:, coord_ind], ydata=teffs, p0=(sigma_0,), bounds=[0.01, 1000])
+
+    sigma = result[0]
+    model = gaussian_1d(coords[:, coord_ind], sigma)
+
+    if plot:
+        plt.scatter(coords[:, coord_ind], teffs)
+        plt.scatter(coords[:, coord_ind], model)
+        plt.show()
+
+    return sigma, amplitude, model
+
+
+def _twoD_Gaussian(coords, sigma_x, sigma_y, amplitude, cutoff=0., offset=0., component=1):
+    """
+    Simple 2D gaussian profile with widths sigma_x and sigma_y and an amplitude.
+    If amplitude is the maximal temperature of a contact component, and the sigmas the width their respective calculated
+    gradients, this function returns the smoothed temperatures of the component.
+
+    Parameters
+    ----------
+    coords: 3D coordinates in roche geometry
+    sigma_x: width of gaussian in x-direction
+    sigma_y: width of gaussian in y-direction
+    amplitude: amplitude of gaussian
+    cutoff: x-displacement of the center of the gaussian (default=0)
+    offset: constant added to the gaussian (default=0)
+    component: which component of the binary (default=1)
+
+    Returns
+    -------
+    2D gaussian profile
+    """
+    y0 = 0.
+    x0 = 0 + cutoff if component == 1 else 1 - cutoff
+    a = 1. / sigma_x ** 2
+    b = 1. / sigma_y ** 2
+    return offset + amplitude * np.exp(- (a * ((coords[:, 0] - x0) ** 2))) * np.exp(- (b * ((coords[:, 1] - y0) ** 2)))
+
+
+def gaussian_smoothing(xyz1, teffs1, xyz2, teffs2, w=0.5, cutoff=0., offset=0.):
+    """
+    Adapts the temperatures of the primary and the secondary to blend nicely using gaussian smoothing.
+    It calculates an average temperature at the neck, appropriate widths of the gaussian profiles (essentially the temp-
+    gradient), and applies the smoothing.
+
+    Parameters
+    ----------
+    xyz1: 3D coordinates in roche geometry of the primary
+    teffs1: teffs of the primary
+    xyz2: 3D coordinates in roche geometry of the secondary
+    teffs2: teffs of the secondary
+    w: contribution of the primary teff vs the secondary (a weighting factor, default=0.5)
+    cutoff: x-displacement of the center of the gaussians toward the center of mass (default=0)
+    offset: constant added to the gaussian (default=0)
+
+    Returns
+    -------
+
+    """
+    coords_neck_1, teffs_neck_1, cond1 = _isolate_neck(xyz1, teffs1, cutoff=cutoff, component=1, plot=False)
+    coords_neck_2, teffs_neck_2, cond2 = _isolate_neck(xyz2, teffs2, cutoff=cutoff, component=2, plot=False)
+
+    x_neck, Tavg = _compute_new_teff_at_neck(coords_neck_1, teffs_neck_1, coords_neck_2, teffs_neck_2, w=w,
+                                             offset=offset)
+
+    sigma_x1 = _compute_sigmax(Tavg, x_neck, x0=0 + cutoff, offset=offset, amplitude=teffs_neck_1.max())
+    sigma_x2 = _compute_sigmax(Tavg, x_neck, x0=1 - cutoff, offset=offset, amplitude=teffs_neck_2.max())
+
+    coords_fit_y1, teffs_fit_y1 = _isolate_sigma_fitting_regions(coords_neck_1, teffs_neck_1, direction='y',
+                                                                 cutoff=cutoff, component=1, plot=False)
+    coords_fit_y2, teffs_fit_y2 = _isolate_sigma_fitting_regions(coords_neck_2, teffs_neck_2, direction='y',
+                                                                 cutoff=cutoff, component=2, plot=False)
+
+    sigma_y1, amplitude_y1, model_y1 = _fit_sigma(coords_fit_y1, teffs_fit_y1, offset, direction='y',
+                                                  component=1, plot=False)
+    sigma_y2, amplitude_y2, model_y2 = _fit_sigma(coords_fit_y2, teffs_fit_y2, offset, direction='y',
+                                                  component=2, plot=False)
+
+    new_teffs1 = _compute_twoD_Gaussian(coords_neck_1, sigma_x1, sigma_y1, teffs_neck_1.max(), cutoff=cutoff,
+                                        offset=offset, component=1)
+    new_teffs2 = _compute_twoD_Gaussian(coords_neck_2, sigma_x2, sigma_y2, teffs_neck_2.max(), cutoff=cutoff,
+                                        offset=offset, component=2)
+
+    # print(cond1, len(cond1), len(new_teffs1))
+    # print(cond2, len(cond2), len(new_teffs2))
+    teffs1[cond1] = new_teffs1
+    teffs2[cond2] = new_teffs2
+
+    return teffs1, teffs2
+
+
+def lateral_transfer(t2s, teffs2, mixing_power, teff_ratio):
+    """
+    Scales the temperatures of the secondary to that of the primary only in a horizontal band the size of the contact's
+    neck. This implies mixing occurs due to mass transfer across the neck.
+    """
+    x2s = t2s[:, 0]
+    y2s = t2s[:, 1]
+    z2s = t2s[:, 2]
+
+    y2s_neck = y2s[x2s < 1]
+    z2s_neck = z2s[x2s < 1]
+    rs_neck = (y2s_neck ** 2 + z2s_neck ** 2) ** 0.5
+    lat = np.min(rs_neck)
+    assert lat == np.min(z2s_neck)
+    filt = (z2s > -lat) & (z2s < lat)  # select band extending the (projected) height of the neck
+    c = (lat - np.abs(z2s[filt])) ** mixing_power
+    latitude_dependence = c / c.max()
+    teffs2[filt] *= 1 + (1 - teff_ratio) * latitude_dependence
+
+    return teffs2
+
+
+def isotropic_transfer(t2s, teffs2, mixing_power, teff_ratio):
+    """
+    Scales the temperatures of the secondary to that of the primary, parametrized by the radial distance from the
+    origin (which is the center of the primary). Implies mixing occurs diffusively from th center of the neck.
+    """
+    d2s = np.sqrt(t2s[:, 0] * t2s[:, 0] + t2s[:, 1] * t2s[:, 1] + t2s[:, 2] * t2s[:, 2])
+    teffs2 *= 1 + (1 - teff_ratio) * (1 - ((d2s - d2s.min()) / (d2s.max() - d2s.min()))) ** mixing_power
+    return teffs2
+
+
+def perfect_transfer(t2s, teff2s, teff_ratio):
+    """
+    Scales the temperatures of the secondary to that of the primary, implying perfect thermal mixing occurred deep in
+    the interior of the stars, and little surface mixing occurs.
+    """
+    teff2s *= 1 / teff_ratio
+    return teff2s
+
+
+def mix_teffs(xyz1, teffs1, xyz2, teffs2, mixing_method='lateral', mixing_power=0.5, teff_ratio=1.):
+    """
+    Applies a temperature mixing, primarily of component 2, according to some mixing method and other parameters.
+    If `mixing_method == 'smoothing'`, simple gaussian smoothing is applied rather than an energy transfer model
+
+    Parameters
+    ----------
+    xyz1: 3D roche coordinates of the primary
+    teffs1: Teffs of the primary (these are only altered by smoothing)
+    xyz2: 3D roche coordinates of the secondary
+    teffs2: Teffs of the secondary (these are altered by all mixing methods)
+    mixing_method: model of energy transfer, 'lateral', 'isotropic', 'spotty', 'perfect' or 'smoothing'
+        (default='lateral')
+    mixing_power: parameter value of the mixing efficiency of the lateral and isotropic mixing model.
+    teff_ratio: ratio of (averaged/global) secondary Teff over (averaged/global) primary Teff (pre-mixing of course).
+        This value determines temperature gradient which drives the mixing. (E.g. if `teff_ratio==1` then no mixing is
+        to be applied, since there's no expected energy transfer)
+
+    Returns
+    -------
+    modified Teffs of the primary, modified Teffs of the secondary
+    """
+    if mixing_method == 'lateral':
+        teffs2 = lateral_transfer(xyz2, teffs2, mixing_power, teff_ratio)
+    elif mixing_method == 'isotropic':
+        teffs2 = isotropic_transfer(xyz2, teffs2, mixing_power, teff_ratio)
+    elif mixing_method == 'perfect':
+        teffs2 = perfect_transfer(xyz2, teffs2, teff_ratio)
+    elif mixing_method == 'smoothing':
+        teffs1, teffs2 = gaussian_smoothing(xyz1, teffs1, xyz2, teffs2)
+    else:
+        raise ValueError('`mixing_method` must be "lateral", "isotropic", "perfect" or "smoothing"')
+    return teffs1, teffs2
diff --git a/phoebe/backend/universe.py b/phoebe/backend/universe.py
index 6424eac15..afee53cad 100644
--- a/phoebe/backend/universe.py
+++ b/phoebe/backend/universe.py
@@ -1581,22 +1581,27 @@ def _fill_teffs(self, mesh=None, ignore_effects=False, **kwargs):
         if mesh is None:
             mesh = self.mesh
 
-        # Now we can compute the local temperatures.
-        # see PHOEBE Legacy scientific reference eq 5.23
-        teffs = self.instantaneous_tpole*mesh.gravs.for_computations**0.25
+        if hasattr(self, 'smoothed_teffs'):
+            teffs = self.smoothed_teffs
+            # print("teffs.median", np.median(teffs))
+        else:
 
-        if not ignore_effects:
-            for feature in self.features:
-                if feature.proto_coords:
+            # Now we can compute the local temperatures.
+            # see PHOEBE Legacy scientific reference eq 5.23
+            teffs = self.instantaneous_tpole*mesh.gravs.for_computations**0.25
 
-                    if self.__class__.__name__ == 'Star_roche_envelope_half' and self.ind_self != self.ind_self_vel:
-                        # then this is the secondary half of a contact envelope
-                        roche_coords_for_computations = np.array([1.0, 0.0, 0.0]) - mesh.roche_coords_for_computations
+            if not ignore_effects:
+                for feature in self.features:
+                    if feature.proto_coords:
+
+                        if self.__class__.__name__ == 'Star_roche_envelope_half' and self.ind_self != self.ind_self_vel:
+                            # then this is the secondary half of a contact envelope
+                            roche_coords_for_computations = np.array([1.0, 0.0, 0.0]) - mesh.roche_coords_for_computations
+                        else:
+                            roche_coords_for_computations = mesh.roche_coords_for_computations
+                        teffs = feature.process_teffs(teffs, roche_coords_for_computations, s=self.polar_direction_xyz, t=self.time)
                     else:
-                        roche_coords_for_computations = mesh.roche_coords_for_computations
-                    teffs = feature.process_teffs(teffs, roche_coords_for_computations, s=self.polar_direction_xyz, t=self.time)
-                else:
-                    teffs = feature.process_teffs(teffs, mesh.coords_for_computations, s=self.polar_direction_xyz, t=self.time)
+                        teffs = feature.process_teffs(teffs, mesh.coords_for_computations, s=self.polar_direction_xyz, t=self.time)
 
         mesh.update_columns(teffs=teffs)
 
diff --git a/phoebe/frontend/default_bundles/default_binary.bundle b/phoebe/frontend/default_bundles/default_binary.bundle
index 6682f0f9f..10c0cfae5 100644
--- a/phoebe/frontend/default_bundles/default_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_binary.bundle
@@ -144,7 +144,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 2.0132751765376384,
+"value": 2.013275176537638,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.437551868254479,
+"value": 4.437551877570185,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185307179586,
+"value": 6.283185,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -460,7 +460,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 0.9988131257959815,
+"value": 0.9988131358058301,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -492,7 +492,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 2.0132751765376384,
+"value": 2.013275176537638,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -557,7 +557,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.437551868254479,
+"value": 4.437551877570185,
 "default_unit": "",
 "limits": [
 null,
@@ -607,7 +607,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185307179586,
+"value": 6.283185,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -808,7 +808,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 0.9988131257959815,
+"value": 0.9988131358058301,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -857,7 +857,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Orbital frequency (sidereal)",
-"value": 6.283185307179586,
+"value": 6.283185,
 "default_unit": "rad / d",
 "limits": [
 null,
@@ -988,7 +988,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Mean anomaly at t0@system",
-"value": 90.0,
+"value": 89.99999559997653,
 "default_unit": "deg",
 "limits": [
 null,
@@ -1245,7 +1245,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@primary@component}",
+"value": "6.283185 / {period@primary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1261,7 +1261,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
+"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1277,7 +1277,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@primary@component}",
+"value": "1.000000 - {irrad_frac_refl_bol@primary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1293,7 +1293,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@secondary@component}",
+"value": "6.283185 / {period@secondary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1309,7 +1309,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
+"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1325,7 +1325,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@secondary@component}",
+"value": "1.000000 - {irrad_frac_refl_bol@secondary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1389,7 +1389,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{period@binary@component} / ((((-1.00000000000000000e+00 * {period@binary@component}) * {dperdt@binary@component}) / 6.28318530717958623e+00) + 1.00000000000000000e+00)",
+"value": "{period@binary@component} / ((((-1.000000 * {period@binary@component}) * {dperdt@binary@component}) / 6.283185307179586231995926937088) + 1.000000000000000000000000000000)",
 "default_unit": "d",
 "constraint_func": "period_anom",
 "constraint_kwargs": {
@@ -1405,7 +1405,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(6.28318530717958623e+00 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
+"value": "(6.283185 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
 "default_unit": "deg",
 "constraint_func": "mean_anom",
 "constraint_kwargs": {
@@ -1463,7 +1463,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@binary@component}",
+"value": "6.283185 / {period@binary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1539,7 +1539,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2.94220621750441933e+03)",
+"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2942.206217504419328179210424423218)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -1561,7 +1561,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{sma@binary@component} / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
+"value": "{sma@binary@component} / ((1.000000 / {q@binary@component}) + 1.000000)",
 "default_unit": "solRad",
 "constraint_func": "comp_sma",
 "constraint_kwargs": {
@@ -1577,7 +1577,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
+"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.000000 / {q@binary@component}) + 1.000000)",
 "default_unit": "solRad",
 "constraint_func": "comp_asini",
 "constraint_kwargs": {
@@ -1657,7 +1657,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)) * 2.94220621750441933e+03)",
+"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.000000 / {q@binary@component}) + 1.000000)) * 2942.206217504419328179210424423218)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -2085,12 +2085,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2291,12 +2290,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2310,12 +2308,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2511,7 +2508,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.15",
+"value": "2.4.16.dev+feature-contacts",
 "copy_for": false,
 "readonly": true,
 "advanced": true,
diff --git a/phoebe/frontend/default_bundles/default_contact_binary.bundle b/phoebe/frontend/default_bundles/default_contact_binary.bundle
index e58fd0f87..5ab782ead 100644
--- a/phoebe/frontend/default_bundles/default_contact_binary.bundle
+++ b/phoebe/frontend/default_bundles/default_contact_binary.bundle
@@ -144,7 +144,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 1.6724563972838287,
+"value": 1.6724563972838384,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -160,7 +160,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (minimum) value of the equivalent radius for the given morphology",
-"value": 1.27254185686813,
+"value": 1.2725418568681297,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.089736153779247,
+"value": 4.089736163094955,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 12.566370614359172,
+"value": 12.56637,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -460,7 +460,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 1.0089067893421308,
+"value": 1.0089067994531355,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -476,7 +476,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Equivalent radius",
-"value": 1.5000000000000122,
+"value": 1.4999999999999996,
 "default_unit": "solRad",
 "limits": [
 1e-06,
@@ -492,7 +492,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 1.672456397283698,
+"value": 1.6724563972838378,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -508,7 +508,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (minimum) value of the equivalent radius for the given morphology",
-"value": 1.27254185686813,
+"value": 1.2725418568681297,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -557,7 +557,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.08973615377924,
+"value": 4.089736163094955,
 "default_unit": "",
 "limits": [
 null,
@@ -607,7 +607,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 12.566370614359172,
+"value": 12.56637,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -808,7 +808,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Mass",
-"value": 1.0089067893421308,
+"value": 1.0089067994531355,
 "default_unit": "solMass",
 "limits": [
 1e-12,
@@ -857,7 +857,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Orbital frequency (sidereal)",
-"value": 12.566370614359172,
+"value": 12.56637,
 "default_unit": "rad / d",
 "limits": [
 null,
@@ -988,7 +988,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "Mean anomaly at t0@system",
-"value": 90.0,
+"value": 89.99999559997653,
 "default_unit": "deg",
 "limits": [
 null,
@@ -1126,12 +1126,56 @@ null
 "Class": "FloatParameter"
 },
 {
+"qualifier": "mixing_enabled",
+"component": "contact_envelope",
+"kind": "envelope",
+"context": "component",
+"description": "Whether to apply thermal mixing of the contact envelope",
+"value": false,
+"copy_for": false,
+"latexfmt": "\\mathrm{{wenabled}}_\\mathrm{{ {component} }}",
+"Class": "BoolParameter"
+},
+{
+"qualifier": "mixing_power",
+"component": "contact_envelope",
+"kind": "envelope",
+"context": "component",
+"description": "Power of the thermal mixing of the envelope",
+"value": 0.5,
+"default_unit": "",
+"limits": [
+null,
+null
+],
+"copy_for": false,
+"latexfmt": "\\mathrm{{w}}_\\mathrm{{ {component} }}",
+"Class": "FloatParameter"
+},
+{
+"qualifier": "mixing_method",
+"component": "contact_envelope",
+"kind": "envelope",
+"context": "component",
+"description": "Method for thermal mixing of the envelope",
+"choices": [
+"lateral",
+"isotropic",
+"perfect",
+"gaussian"
+],
+"value": "lateral",
+"copy_for": false,
+"latexfmt": "\\mathrm{{smethod}}_\\mathrm{{ {component} }}",
+"Class": "ChoiceParameter"
+},
+{
 "qualifier": "fillout_factor",
 "component": "contact_envelope",
 "kind": "envelope",
 "context": "component",
 "description": "Fillout-factor of the envelope",
-"value": 0.6417897080770861,
+"value": 0.6417897080770943,
 "default_unit": "",
 "limits": [
 0.0,
@@ -1147,7 +1191,7 @@ null
 "kind": "envelope",
 "context": "component",
 "description": "Potential of the envelope (from the primary component's reference)",
-"value": 3.4013774072298815,
+"value": 3.401377407229877,
 "default_unit": "",
 "limits": [
 0.0,
@@ -1211,7 +1255,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "ratio between equivalent radii of children stars",
-"value": 1.0000000000000082,
+"value": 0.9999999999999997,
 "default_unit": "",
 "limits": [
 0.0,
@@ -1227,7 +1271,7 @@ null
 "kind": "orbit",
 "context": "component",
 "description": "sum of fractional equivalent radii of children stars",
-"value": 0.8955223880597052,
+"value": 0.8955223880597013,
 "default_unit": "",
 "limits": [
 0.0,
@@ -1324,7 +1368,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@primary@component}",
+"value": "6.283185 / {period@primary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1340,7 +1384,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
+"value": "log10((({mass@primary@component} / ({requiv@primary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1356,7 +1400,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@primary@component}",
+"value": "1.000000 - {irrad_frac_refl_bol@primary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1372,7 +1416,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@secondary@component}",
+"value": "6.283185 / {period@secondary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1388,7 +1432,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
+"value": "log10((({mass@secondary@component} / ({requiv@secondary@component} ** 2.000000)) * 2942.206218) * 9.319541)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -1404,7 +1448,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@secondary@component}",
+"value": "1.000000 - {irrad_frac_refl_bol@secondary@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -1468,7 +1512,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{period@binary@component} / ((((-1.00000000000000000e+00 * {period@binary@component}) * {dperdt@binary@component}) / 6.28318530717958623e+00) + 1.00000000000000000e+00)",
+"value": "{period@binary@component} / ((((-1.000000 * {period@binary@component}) * {dperdt@binary@component}) / 6.283185307179586231995926937088) + 1.000000000000000000000000000000)",
 "default_unit": "d",
 "constraint_func": "period_anom",
 "constraint_kwargs": {
@@ -1484,7 +1528,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(6.28318530717958623e+00 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
+"value": "(6.283185 * ({t0@system} - {t0_perpass@binary@component})) / {period@binary@component}",
 "default_unit": "deg",
 "constraint_func": "mean_anom",
 "constraint_kwargs": {
@@ -1542,7 +1586,7 @@ null
 "kind": "orbit",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@binary@component}",
+"value": "6.283185 / {period@binary@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -1698,7 +1742,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2.94220621750441933e+03)",
+"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ({q@binary@component} + 1.000000)) * 2942.206217504419328179210424423218)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -1720,7 +1764,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "{sma@binary@component} / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
+"value": "{sma@binary@component} / ((1.000000 / {q@binary@component}) + 1.000000)",
 "default_unit": "solRad",
 "constraint_func": "comp_sma",
 "constraint_kwargs": {
@@ -1736,7 +1780,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)",
+"value": "({sma@binary@component} * (sin({incl@binary@component}))) / ((1.000000 / {q@binary@component}) + 1.000000)",
 "default_unit": "solRad",
 "constraint_func": "comp_asini",
 "constraint_kwargs": {
@@ -1832,7 +1876,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "(3.94784176043574320e+01 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.00000000000000000e+00 / {q@binary@component}) + 1.00000000000000000e+00)) * 2.94220621750441933e+03)",
+"value": "(39.478418 * ({sma@binary@component} ** 3.000000)) / ((({period@binary@component} ** 2.000000) * ((1.000000 / {q@binary@component}) + 1.000000)) * 2942.206217504419328179210424423218)",
 "default_unit": "solMass",
 "constraint_func": "mass",
 "constraint_kwargs": {
@@ -2276,12 +2320,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -2513,12 +2556,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2532,12 +2574,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -2743,7 +2784,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.15",
+"value": "2.4.16.dev+feature-contacts",
 "copy_for": false,
 "readonly": true,
 "advanced": true,
diff --git a/phoebe/frontend/default_bundles/default_star.bundle b/phoebe/frontend/default_bundles/default_star.bundle
index a93ac8444..f85293fec 100644
--- a/phoebe/frontend/default_bundles/default_star.bundle
+++ b/phoebe/frontend/default_bundles/default_star.bundle
@@ -144,7 +144,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Critical (maximum) value of the equivalent radius for the given morphology",
-"value": 3.429265859647371,
+"value": 3.4292622920441116,
 "default_unit": "solRad",
 "limits": [
 0.0,
@@ -209,7 +209,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "logg at requiv",
-"value": 4.438067627303133,
+"value": 4.438067632266453,
 "default_unit": "",
 "limits": [
 null,
@@ -259,7 +259,7 @@ null
 "kind": "star",
 "context": "component",
 "description": "Rotation frequency (wrt the sky)",
-"value": 6.283185307179586,
+"value": 6.283185,
 "default_unit": "rad / d",
 "limits": [
 0.0,
@@ -489,7 +489,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "6.28318530717958623e+00 / {period@starA@component}",
+"value": "6.283185 / {period@starA@component}",
 "default_unit": "rad / d",
 "constraint_func": "freq",
 "constraint_kwargs": {
@@ -505,7 +505,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "log10((({mass@starA@component} / ({requiv@starA@component} ** 2.000000)) * 2.94220621750441933e+03) * 9.31954089506172778e+00)",
+"value": "log10((({mass@starA@component} / ({requiv@starA@component} ** 2.000000)) * 2942.206218) * 9.319541)",
 "default_unit": "",
 "constraint_func": "logg",
 "constraint_kwargs": {
@@ -521,7 +521,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "1.00000000000000000e+00 - {irrad_frac_refl_bol@starA@component}",
+"value": "1.000000 - {irrad_frac_refl_bol@starA@component}",
 "default_unit": "",
 "constraint_func": "irrad_frac",
 "constraint_kwargs": {
@@ -537,7 +537,7 @@ null
 "kind": "star",
 "context": "constraint",
 "description": "expression that determines the constraint",
-"value": "8.14885676770049971e-01 * (((2.94220621750441933e+03 * {mass@starA@component}) * (({period@starA@component} / 6.283185) ** 2.00000000000000000e+00)) ** 3.33333333333333315e-01)",
+"value": "0.814886 * (((2942.206218 * {mass@starA@component}) * (({period@starA@component} / 6.283185) ** 2.000000)) ** 0.333333)",
 "default_unit": "solRad",
 "constraint_func": "requiv_single_max",
 "constraint_kwargs": {
@@ -863,12 +863,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": {
@@ -1020,12 +1019,11 @@ null
 "context": "compute",
 "description": "Atmosphere table",
 "choices": [
-"phoenix",
 "extern_atmx",
-"extern_planckint",
-"blackbody",
 "phoenix",
-"ck2004"
+"extern_planckint",
+"ck2004",
+"blackbody"
 ],
 "value": "ck2004",
 "copy_for": false,
@@ -1136,7 +1134,7 @@ null
 "qualifier": "phoebe_version",
 "context": "setting",
 "description": "Version of PHOEBE",
-"value": "2.4.15",
+"value": "2.4.16.dev+feature-contacts",
 "copy_for": false,
 "readonly": true,
 "advanced": true,
diff --git a/phoebe/parameters/component.py b/phoebe/parameters/component.py
index 257208974..563566e9f 100644
--- a/phoebe/parameters/component.py
+++ b/phoebe/parameters/component.py
@@ -275,10 +275,13 @@ def envelope(component, **kwargs):
     * `abun` (float, optional): abundance/metallicity.
     * `fillout_factor` (float, optional): fillout-factor of the envelope.
     * `pot` (float, optional): potential of the envelope.
-    * `pot_min` (float, optional): critical (minimum) value of the potential to
-        remain a contact.
-    * `pot_max` (float, optional): critical (maximum) value of the potential to
-        remain a contact.
+    * `pot_min` (float, optional): critical (minimum) value of the potential to remain a contact.
+    * `pot_max` (float, optional): critical (maximum) value of the potential to remain a contact.
+    * `mixing_enabled` (float, optional): whether to allow thermal mixing of the surface (default = False)
+    * `mixing_method` (str, optional): method to use to perform the thermal mixing, choose between [lateral, isotropic,
+            perfect, spotty] (default = 'lateral')
+    * `mixing_power` (float, optional): thermal mixing strength, free parameter of the model. Generally, higher values
+            mean more efficient mixing (default = 0.5)
 
     Returns
     --------
@@ -298,6 +301,9 @@ def envelope(component, **kwargs):
     # params += [FloatParameter(qualifier='frac_lost_bol', value=kwargs.get('frac_lost_bol', 1.0), default_unit=u.dimensionless_unscaled, limits=(0.0, 1.0), description='ratio of incident bolometric light that is lost/ignored')]
 
 
+    params += [BoolParameter(qualifier='mixing_enabled', latexfmt=r'\mathrm{{wenabled}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_enabled', False), description='Whether to apply thermal mixing of the contact envelope')]
+    params += [FloatParameter(qualifier='mixing_power', latexfmt=r'\mathrm{{w}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_power', 0.5), default_unit=u.dimensionless_unscaled, limits=(None,None), description='Power of the thermal mixing of the envelope')]
+    params += [ChoiceParameter(qualifier='mixing_method', latexfmt=r'\mathrm{{smethod}}_\mathrm{{ {component} }}', value=kwargs.get('mixing_method', 'lateral'), choices=['lateral', 'isotropic', 'perfect', 'gaussian'], description='Method for thermal mixing of the envelope')]
     params += [FloatParameter(qualifier='fillout_factor', latexfmt=r'\mathrm{{FF}}_\mathrm{{ {component} }}', value=kwargs.get('fillout_factor', 0.5), default_unit=u.dimensionless_unscaled, limits=(0.0,1.0), description='Fillout-factor of the envelope')]
     params += [FloatParameter(qualifier='pot', latexfmt=r'\Omega_\mathrm{{ {component} }}', value=kwargs.get('pot', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Potential of the envelope (from the primary component\'s reference)')]
     params += [FloatParameter(qualifier='pot_min', latexfmt=r'\Omega_\mathrm{{ min,  {component} }}', value=kwargs.get('pot_min', 3.5), default_unit=u.dimensionless_unscaled, limits=(0.0,None), description='Critical (minimum) value of the potential to remain a contact')]
diff --git a/phoebe/parameters/parameters.py b/phoebe/parameters/parameters.py
index 128bdbd85..3f3f66e59 100644
--- a/phoebe/parameters/parameters.py
+++ b/phoebe/parameters/parameters.py
@@ -208,7 +208,8 @@
                       'mass', 'dpdt', 'per0',
                       'dperdt', 'ecc', 'deccdt', 't0_perpass', 't0_supconj',
                       't0_ref', 'mean_anom', 'q', 'sma', 'asini', 'ecosw', 'esinw',
-                      'teffratio', 'requivratio', 'requivsumfrac'
+                      'teffratio', 'requivratio', 'requivsumfrac',
+                      "mixing_enabled", "mixing_power", "mixing_method"
                       ]
 
 # from dataset:
diff --git a/tests/tests/test_contacts/test_mixing.py b/tests/tests/test_contacts/test_mixing.py
new file mode 100644
index 000000000..9769e21d3
--- /dev/null
+++ b/tests/tests/test_contacts/test_mixing.py
@@ -0,0 +1,23 @@
+import phoebe
+import numpy as np
+
+
+def test_mixing(verbose=False):
+    b = phoebe.default_contact_binary()
+    b.add_dataset('mesh', compute_times=[0])
+    b.set_value('teff@secondary', 3000)  # so we have different temperatures
+    b.set_value('mixing_enabled', True)
+    b.set_value('mixing_method', 'perfect')  # since other mixing methods are arbitrarily parametrized, test this case only
+    b.set_value('columns', '*')  # to populate 'teffs'
+
+    assert b.run_checks().passed
+
+    b.run_compute()
+
+    assert (res := abs(1 - np.mean(b['teffs@secondary'].value) / np.mean(b['teffs@primary'].value))) <= 0.01
+    if verbose:
+        print('relative teff difference', res)
+
+
+if __name__ == '__main__':
+    test_mixing(verbose=True)