diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py
index 2627fdc698fff114fc42d0bf9d89914551b833e4..cc903044b59b2a33f04d50292c0d8aabfdcbf486 100644
--- a/demos/nonlinear_critical_filter.py
+++ b/demos/nonlinear_critical_filter.py
@@ -59,9 +59,8 @@ if __name__ == "__main__":
 
     d_space = MeasurementOperator.target
 
-    Projection = ift.PowerProjectionOperator(domain=h_space,
-                                             power_space=p_space)
-    power = Projection.adjoint_times(ift.exp(0.5*log_p))
+    Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
+    power = Distributor(ift.exp(0.5*log_p))
     # Creating the mock data
     true_sky = nonlinearity(HT(power*sh))
     noiseless_data = MeasurementOperator(true_sky)
@@ -75,7 +74,7 @@ if __name__ == "__main__":
 
     m0 = ift.power_synthesize(ift.Field(p_space, val=1e-7))
     t0 = ift.Field(p_space, val=-4.)
-    power0 = Projection.adjoint_times(ift.exp(0.5 * t0))
+    power0 = Distributor.times(ift.exp(0.5 * t0))
 
     IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
                                      tol_abs_gradnorm=1e-3)
@@ -87,7 +86,7 @@ if __name__ == "__main__":
     inverter = ift.ConjugateGradient(controller=ICI)
 
     for i in range(20):
-        power0 = Projection.adjoint_times(ift.exp(0.5*t0))
+        power0 = Distributor(ift.exp(0.5*t0))
         map0_energy = ift.library.NonlinearWienerFilterEnergy(
             m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
             inverter=inverter)
@@ -103,7 +102,7 @@ if __name__ == "__main__":
         power0_energy = ift.library.NonlinearPowerEnergy(
             position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
             Instrument=MeasurementOperator, nonlinearity=nonlinearity,
-            Projection=Projection, sigma=1., samples=2, inverter=inverter)
+            Distributor=Distributor, sigma=1., samples=2, inverter=inverter)
 
         power0_energy = minimizer(power0_energy)[0]
 
diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py
index 8c6e11105d556ab1b75bcfedc04b67718703a7ba..8dd1fb854e87ae535d43b116a838a9311fad5753 100644
--- a/demos/wiener_filter_via_hamiltonian.py
+++ b/demos/wiener_filter_via_hamiltonian.py
@@ -70,7 +70,8 @@ if __name__ == "__main__":
     # Generate plots
     zmax = max(ht(sh).max(), ht(m).max())
     zmin = min(ht(sh).min(), ht(m).min())
-    plotdict = {"zmax": zmax, "zmin": zmin, "colormap": "Planck-like"}
+    plotdict = {"zmax": zmax, "zmin": zmin,"colormap": "Planck-like"}
+    plotdict2 = {"colormap": "Planck-like"}
     ift.plot(ht(sh), name="mock_signal.png", **plotdict)
     ift.plot(ht(m), name="reconstruction.png", **plotdict)
 
@@ -79,10 +80,10 @@ if __name__ == "__main__":
     sample_mean = ift.Field.zeros(s_space)
 
     mean, variance = ift.probe_with_posterior_samples(curv, ht, 50)
-    ift.plot(variance, name="posterior_variance.png", **plotdict)
+    ift.plot(variance, name="posterior_variance.png", **plotdict2)
     ift.plot(mean+ht(m), name="posterior_mean.png", **plotdict)
 
     # try to do the same with diagonal probing
     variance = ift.probe_diagonal(ht*curv.inverse*ht.adjoint, 100)
-    #sm = ift.FFTSmoothingOperator(s_space, sigma=0.015)
-    ift.plot(variance, name="posterior_variance2.png", **plotdict)
+    #sm = ift.FFTSmoothingOperator(s_space, sigma=0.005)
+    ift.plot(variance, name="posterior_variance2.png", **plotdict2)
diff --git a/docs/source/concepts/operators.rst b/docs/source/concepts/operators.rst
index f6ebfcf6c08611039c36ab1ca1180861b6379b4d..21a95926665e3127576298820e260bd802c1f3e1 100644
--- a/docs/source/concepts/operators.rst
+++ b/docs/source/concepts/operators.rst
@@ -12,7 +12,7 @@ Description of Operators
     ScalingOperator                     <../mod/nifty4.operators.scaling_operator>
     DiagonalOperator                    <../mod/nifty4.operators.diagonal_operator>
     HarmonicTransformOperator           <../mod/nifty4.operators.harmonic_transform_operator>
-    PowerProjectionOperator             <../mod/nifty4.operators.power_projection_operator>
+    PowerDistributor                    <../mod/nifty4.operators.power_distributor>
 
 .. toctree::
     :maxdepth: 1
@@ -26,7 +26,7 @@ Description of Operators
 
 ..    Adjoint                             <../mod/nifty4.operators.adjoint_operator>
 ..    Chain                               <../mod/nifty4.operators.chain_operator>
-..    DOF Projection                      <../mod/nifty4.operators.dof_projection_operator>
+..    DOF Distributor                     <../mod/nifty4.operators.dof_distributor>
 ..    Inverse                             <../mod/nifty4.operators.inverse_operator>
 ..    Laplace                             <../mod/nifty4.operators.linear_operator>
 ..    Sum                                 <../mod/nifty4.operators.sum_operator>
diff --git a/nifty4/library/critical_power_energy.py b/nifty4/library/critical_power_energy.py
index bc08071314556928f6baa177aa419e3d9f46473b..f94ceb16e8dc9d8df8fc786e3f311f76b83892c5 100644
--- a/nifty4/library/critical_power_energy.py
+++ b/nifty4/library/critical_power_energy.py
@@ -18,7 +18,7 @@
 
 from ..minimization.energy import Energy
 from ..operators.smoothness_operator import SmoothnessOperator
-from ..operators.power_projection_operator import PowerProjectionOperator
+from ..operators.power_distributor import PowerDistributor
 from .critical_power_curvature import CriticalPowerCurvature
 from ..utilities import memo
 from .. import Field, exp
@@ -86,17 +86,17 @@ class CriticalPowerEnergy(Energy):
         self._inverter = inverter
 
         if w is None:
-            P = PowerProjectionOperator(domain=self.m.domain,
-                                        power_space=self.position.domain[0])
+            Dist = PowerDistributor(target=self.m.domain,
+                                    power_space=self.position.domain[0])
             if self.D is not None:
                 w = Field.zeros(self.position.domain, dtype=self.m.dtype)
                 for i in range(self.samples):
                     sample = self.D.draw_sample() + self.m
-                    w += P(abs(sample)**2)
+                    w += Dist.adjoint_times(abs(sample)**2)
 
                 w *= 1./self.samples
             else:
-                w = P(abs(self.m)**2)
+                w = Dist.adjoint_times(abs(self.m)**2)
         self._w = w
 
         self._theta = exp(-self.position) * (self.q + self._w*0.5)
diff --git a/nifty4/library/noise_energy.py b/nifty4/library/noise_energy.py
index 465b3a621d56cae0085ae1d028dff33dd618b874..419df97f3d3229d9161bb38ffc0709e2b568935d 100644
--- a/nifty4/library/noise_energy.py
+++ b/nifty4/library/noise_energy.py
@@ -25,7 +25,7 @@ from ..minimization.energy import Energy
 
 class NoiseEnergy(Energy):
     def __init__(self, position, d, xi, D, t, ht, Instrument,
-                 nonlinearity, alpha, q, Projection, samples=3,
+                 nonlinearity, alpha, q, Distribution, samples=3,
                  xi_sample_list=None, inverter=None):
         super(NoiseEnergy, self).__init__(position=position)
         self.xi = xi
@@ -40,8 +40,8 @@ class NoiseEnergy(Energy):
 
         self.alpha = alpha
         self.q = q
-        self.Projection = Projection
-        self.power = self.Projection.adjoint_times(exp(0.5 * self.t))
+        self.Distribution = Distribution
+        self.power = self.Distribution(exp(0.5 * self.t))
         if xi_sample_list is None:
             if samples is None or samples == 0:
                 xi_sample_list = [xi]
@@ -51,7 +51,7 @@ class NoiseEnergy(Energy):
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
 
-        A = Projection.adjoint_times(exp(.5*self.t))
+        A = Distribution(exp(.5*self.t))
 
         self._gradient = None
         for sample in self.xi_sample_list:
@@ -81,7 +81,7 @@ class NoiseEnergy(Energy):
         return self.__class__(
             position, self.d, self.xi, self.D, self.t, self.ht,
             self.Instrument, self.nonlinearity, self.alpha, self.q,
-            self.Projection, xi_sample_list=self.xi_sample_list,
+            self.Distribution, xi_sample_list=self.xi_sample_list,
             samples=self.samples, inverter=self.inverter)
 
     @property
diff --git a/nifty4/library/nonlinear_power_energy.py b/nifty4/library/nonlinear_power_energy.py
index 670757d1f64c25efa18dc471e02eafe41c3b2034..dc78a904b8f398c521d08f8beb2513be065142b0 100644
--- a/nifty4/library/nonlinear_power_energy.py
+++ b/nifty4/library/nonlinear_power_energy.py
@@ -52,7 +52,7 @@ class NonlinearPowerEnergy(Energy):
     """
 
     def __init__(self, position, d, N, xi, D, ht, Instrument, nonlinearity,
-                 Projection, sigma=0., samples=3, xi_sample_list=None,
+                 Distribution, sigma=0., samples=3, xi_sample_list=None,
                  inverter=None):
         super(NonlinearPowerEnergy, self).__init__(position)
         self.xi = xi
@@ -64,7 +64,7 @@ class NonlinearPowerEnergy(Energy):
         self.ht = ht
         self.Instrument = Instrument
         self.nonlinearity = nonlinearity
-        self.Projection = Projection
+        self.Distribution = Distribution
         self.sigma = sigma
         if xi_sample_list is None:
             if samples is None or samples == 0:
@@ -75,7 +75,7 @@ class NonlinearPowerEnergy(Energy):
         self.xi_sample_list = xi_sample_list
         self.inverter = inverter
 
-        A = Projection.adjoint_times(exp(.5 * position))
+        A = Distribution(exp(.5 * position))
         map_s = self.ht(A * xi)
         Tpos = self.T(position)
 
@@ -83,7 +83,7 @@ class NonlinearPowerEnergy(Energy):
         for xi_sample in self.xi_sample_list:
             map_s = self.ht(A * xi_sample)
             LinR = LinearizedPowerResponse(
-                self.Instrument, self.nonlinearity, self.ht, self.Projection,
+                self.Instrument, self.nonlinearity, self.ht, self.Distribution,
                 self.position, xi_sample)
 
             residual = self.d - \
@@ -106,7 +106,7 @@ class NonlinearPowerEnergy(Energy):
     def at(self, position):
         return self.__class__(position, self.d, self.N, self.xi, self.D,
                               self.ht, self.Instrument, self.nonlinearity,
-                              self.Projection, sigma=self.sigma,
+                              self.Distribution, sigma=self.sigma,
                               samples=len(self.xi_sample_list),
                               xi_sample_list=self.xi_sample_list,
                               inverter=self.inverter)
@@ -124,5 +124,5 @@ class NonlinearPowerEnergy(Energy):
     def curvature(self):
         return NonlinearPowerCurvature(
             self.position, self.ht, self.Instrument, self.nonlinearity,
-            self.Projection, self.N, self.T, self.xi_sample_list,
+            self.Distribution, self.N, self.T, self.xi_sample_list,
             self.inverter)
diff --git a/nifty4/library/response_operators.py b/nifty4/library/response_operators.py
index ec7edfe783b18957ec34208fe88c6f9c45f0f3f9..b1fbec30767c277751ecad31e5feedb77b3424fd 100644
--- a/nifty4/library/response_operators.py
+++ b/nifty4/library/response_operators.py
@@ -23,8 +23,9 @@ def LinearizedSignalResponse(Instrument, nonlinearity, ht, power, m):
     return Instrument * nonlinearity.derivative(m) * ht * power
 
 
-def LinearizedPowerResponse(Instrument, nonlinearity, ht, Projection, tau, xi):
+def LinearizedPowerResponse(Instrument, nonlinearity, ht, Distribution, tau,
+        xi):
     power = exp(0.5*tau)
-    position = ht(Projection.adjoint_times(power)*xi)
+    position = ht(Distribution(power)*xi)
     linearization = nonlinearity.derivative(position)
-    return 0.5*Instrument*linearization*ht*xi*Projection.adjoint*power
+    return 0.5*Instrument*linearization*ht*xi*Distribution*power
diff --git a/nifty4/operators/__init__.py b/nifty4/operators/__init__.py
index 200776217b4501a340be4404b11efe1182ac62fc..d2d1207802ea3751ad18054b8dc0cb1583e6072b 100644
--- a/nifty4/operators/__init__.py
+++ b/nifty4/operators/__init__.py
@@ -7,10 +7,10 @@ from .fft_operator import FFTOperator
 from .fft_smoothing_operator import FFTSmoothingOperator
 from .geometry_remover import GeometryRemover
 from .laplace_operator import LaplaceOperator
-from .power_projection_operator import PowerProjectionOperator
+from .power_distributor import PowerDistributor
 from .inversion_enabler import InversionEnabler
 
 __all__ = ["LinearOperator", "EndomorphicOperator", "ScalingOperator",
            "DiagonalOperator", "HarmonicTransformOperator", "FFTOperator",
            "FFTSmoothingOperator", "GeometryRemover",
-           "LaplaceOperator", "PowerProjectionOperator", "InversionEnabler"]
+           "LaplaceOperator", "PowerDistributor", "InversionEnabler"]
diff --git a/nifty4/operators/dof_projection_operator.py b/nifty4/operators/dof_distributor.py
similarity index 79%
rename from nifty4/operators/dof_projection_operator.py
rename to nifty4/operators/dof_distributor.py
index 1ab10f26d7a63d3aa6fab23cf92d863a74b72872..ebd8b8277449473203348a44e9f940dce164ca36 100644
--- a/nifty4/operators/dof_projection_operator.py
+++ b/nifty4/operators/dof_distributor.py
@@ -25,15 +25,15 @@ from .. import dobj
 from ..domains.dof_space import DOFSpace
 
 
-class DOFProjectionOperator(LinearOperator):
-    def __init__(self, dofdex, domain=None, space=None):
-        super(DOFProjectionOperator, self).__init__()
+class DOFDistributor(LinearOperator):
+    def __init__(self, dofdex, target=None, space=None):
+        super(DOFDistributor, self).__init__()
 
-        if domain is None:
-            domain = dofdex.domain
-        self._domain = DomainTuple.make(domain)
-        space = infer_space(self._domain, space)
-        partner = self._domain[space]
+        if target is None:
+            target = dofdex.domain
+        self._target = DomainTuple.make(target)
+        space = infer_space(self._target, space)
+        partner = self._target[space]
         if not isinstance(dofdex, Field):
             raise TypeError("dofdex must be a Field")
         if not len(dofdex.domain) == 1:
@@ -64,41 +64,41 @@ class DOFProjectionOperator(LinearOperator):
 
     def _init2(self, dofdex, space, other_space):
         self._space = space
-        tgt = list(self._domain)
-        tgt[self._space] = other_space
-        self._target = DomainTuple.make(tgt)
+        dom = list(self._target)
+        dom[self._space] = other_space
+        self._domain = DomainTuple.make(dom)
 
         if dobj.default_distaxis() in self._domain.axes[self._space]:
             dofdex = dobj.local_data(dofdex)
         else:  # dofdex must be available fully on every task
             dofdex = dobj.to_global_data(dofdex)
         self._dofdex = dofdex.ravel()
-        firstaxis = self._domain.axes[self._space][0]
-        lastaxis = self._domain.axes[self._space][-1]
-        arrshape = dobj.local_shape(self._domain.shape, 0)
+        firstaxis = self._target.axes[self._space][0]
+        lastaxis = self._target.axes[self._space][-1]
+        arrshape = dobj.local_shape(self._target.shape, 0)
         presize = np.prod(arrshape[0:firstaxis], dtype=np.int)
         postsize = np.prod(arrshape[lastaxis+1:], dtype=np.int)
-        self._hshape = (presize, self._target[self._space].shape[0], postsize)
+        self._hshape = (presize, self._domain[self._space].shape[0], postsize)
         self._pshape = (presize, self._dofdex.size, postsize)
 
-    def _times(self, x):
+    def _adjoint_times(self, x):
         arr = dobj.local_data(x.val)
         arr = arr.reshape(self._pshape)
         oarr = np.zeros(self._hshape, dtype=x.dtype)
         np.add.at(oarr, (slice(None), self._dofdex, slice(None)), arr)
         if dobj.distaxis(x.val) in x.domain.axes[self._space]:
-            oarr = dobj.np_allreduce_sum(oarr).reshape(self._target.shape)
-            res = Field(self._target, dobj.from_global_data(oarr))
+            oarr = dobj.np_allreduce_sum(oarr).reshape(self._domain.shape)
+            res = Field(self._domain, dobj.from_global_data(oarr))
         else:
-            oarr = oarr.reshape(dobj.local_shape(self._target.shape,
+            oarr = oarr.reshape(dobj.local_shape(self._domain.shape,
                                                  dobj.distaxis(x.val)))
-            res = Field(self._target,
-                        dobj.from_local_data(self._target.shape, oarr,
+            res = Field(self._domain,
+                        dobj.from_local_data(self._domain.shape, oarr,
                                              dobj.default_distaxis()))
         return res
 
-    def _adjoint_times(self, x):
-        res = Field.empty(self._domain, dtype=x.dtype)
+    def _times(self, x):
+        res = Field.empty(self._target, dtype=x.dtype)
         if dobj.distaxis(x.val) in x.domain.axes[self._space]:
             arr = dobj.to_global_data(x.val)
         else:
diff --git a/nifty4/operators/power_projection_operator.py b/nifty4/operators/power_distributor.py
similarity index 78%
rename from nifty4/operators/power_projection_operator.py
rename to nifty4/operators/power_distributor.py
index c6a01aeee1c1c7919524c38cf91a4e58a63bd5ee..5810e2e0d38f79ddd3a5cf0eb261f60deba98b25 100644
--- a/nifty4/operators/power_projection_operator.py
+++ b/nifty4/operators/power_distributor.py
@@ -16,20 +16,20 @@
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
-from .dof_projection_operator import DOFProjectionOperator
+from .dof_distributor import DOFDistributor
 from ..domain_tuple import DomainTuple
 from ..utilities import infer_space
 from ..domains.power_space import PowerSpace
 
 
-class PowerProjectionOperator(DOFProjectionOperator):
-    def __init__(self, domain, power_space=None, space=None):
+class PowerDistributor(DOFDistributor):
+    def __init__(self, target, power_space=None, space=None):
         # Initialize domain and target
-        self._domain = DomainTuple.make(domain)
-        self._space = infer_space(self._domain, space)
-        hspace = self._domain[self._space]
+        self._target = DomainTuple.make(target)
+        self._space = infer_space(self._target, space)
+        hspace = self._target[self._space]
         if not hspace.harmonic:
-            raise ValueError("Operator acts on harmonic spaces only")
+            raise ValueError("Operator requires harmonic target space")
         if power_space is None:
             power_space = PowerSpace(hspace)
         else:
diff --git a/nifty4/sugar.py b/nifty4/sugar.py
index f3204c3ac63e8763d03bfc9ab0668d07df98e5d2..16f2dd70537ea997a18323818aa2c6ccb2667fb3 100644
--- a/nifty4/sugar.py
+++ b/nifty4/sugar.py
@@ -21,7 +21,7 @@ from .domains.structured_domain import StructuredDomain
 from .domains.power_space import PowerSpace
 from .field import Field, sqrt
 from .operators.diagonal_operator import DiagonalOperator
-from .operators.power_projection_operator import PowerProjectionOperator
+from .operators.power_distributor import PowerDistributor
 from .operators.harmonic_transform_operator import HarmonicTransformOperator
 from .domain_tuple import DomainTuple
 from . import dobj, utilities
@@ -45,8 +45,8 @@ def PS_field(pspace, func, dtype=None):
 
 def _single_power_analyze(field, idx, binbounds):
     power_domain = PowerSpace(field.domain[idx], binbounds)
-    ppo = PowerProjectionOperator(field.domain, power_domain, idx)
-    return ppo(field.weight(1)).weight(-1)  # divides by bin size
+    pd = PowerDistributor(field.domain, power_domain, idx)
+    return pd.adjoint_times(field.weight(1)).weight(-1)  # divides by bin size
 
 
 def power_analyze(field, spaces=None, binbounds=None,
@@ -125,8 +125,8 @@ def power_synthesize_nonrandom(field, spaces=None):
     spec = sqrt(field)
     for i in spaces:
         result_domain[i] = field.domain[i].harmonic_partner
-        ppo = PowerProjectionOperator(result_domain, field.domain[i], i)
-        spec = ppo.adjoint_times(spec)
+        pd = PowerDistributor(result_domain, field.domain[i], i)
+        spec = pd(spec)
 
     return spec
 
@@ -203,7 +203,7 @@ def create_power_field(domain, power_spectrum, dtype=None):
         power_domain = PowerSpace(domain)
         fp = PS_field(power_domain, power_spectrum, dtype)
 
-    return PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
+    return PowerDistributor(domain, power_domain)(fp)
 
 
 def create_power_operator(domain, power_spectrum, space=None, dtype=None):
diff --git a/test/test_energies/test_map.py b/test/test_energies/test_map.py
index e00bda8b550bd80a2149b8af14530e84b562a469..21f296aaa0b10ce1cf34c877c062561b8ea86eea 100644
--- a/test/test_energies/test_map.py
+++ b/test/test_energies/test_map.py
@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2017 Max-Planck-Society
+# Copyright(C) 2013-2018 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
@@ -39,12 +39,12 @@ class Energy_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
         s0 = xi0 * A
         Instrument = ift.ScalingOperator(10., space)
@@ -85,12 +85,12 @@ class Energy_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
         s = ht(xi0 * A)
         R = ift.ScalingOperator(10., space)
@@ -129,12 +129,12 @@ class Curvature_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
         s0 = xi0 * A
         Instrument = ift.ScalingOperator(10., space)
@@ -178,12 +178,12 @@ class Curvature_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal')
         s = ht(xi0 * A)
         R = ift.ScalingOperator(10., space)
diff --git a/test/test_energies/test_noise.py b/test/test_energies/test_noise.py
index 30d6d256fd64752be1ed4d2bf2c2cfde9a414329..dde5bf9006e318664f55e21db1ce92f78ef6a106 100644
--- a/test/test_energies/test_noise.py
+++ b/test/test_energies/test_noise.py
@@ -11,7 +11,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright(C) 2013-2017 Max-Planck-Society
+# Copyright(C) 2013-2018 Max-Planck-Society
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
 # and financially supported by the Studienstiftung des deutschen Volkes.
@@ -41,12 +41,12 @@ class Noise_Energy_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, target=space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=False)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         tau = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(tau))
+        A = Dist(ift.sqrt(tau))
         var = 1.
         n = ift.Field.from_random(
             domain=space,
@@ -86,7 +86,7 @@ class Noise_Energy_Tests(unittest.TestCase):
 
         energy0 = ift.library.NoiseEnergy(
             position=eta0, d=d, xi=xi, D=D, t=tau, Instrument=R,
-            alpha=alpha, q=q, Projection=P, nonlinearity=f,
+            alpha=alpha, q=q, Distribution=Dist, nonlinearity=f,
             ht=ht, samples=10)
         energy1 = energy0.at(eta1)
 
diff --git a/test/test_energies/test_power.py b/test/test_energies/test_power.py
index 2386232602595a3b9de759d7eaf8360d30528ef3..31a093de4d3455470330ca7fa21dfb656871f40b 100644
--- a/test/test_energies/test_power.py
+++ b/test/test_energies/test_power.py
@@ -39,13 +39,13 @@ class Energy_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=True)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 64 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
         tau0 = ift.log(pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal', std=.01)
         N = ift.DiagonalOperator(n**2)
         s = xi * A
@@ -91,12 +91,12 @@ class Energy_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=True)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 1 / (1 + k**2)**dim
         tau0 = ift.PS_field(pspace, pspec)
-        A = P.adjoint_times(ift.sqrt(tau0))
+        A = Dist(ift.sqrt(tau0))
         n = ift.Field.from_random(domain=space, random_type='normal')
         s = ht(xi * A)
         R = ift.ScalingOperator(10., space)
@@ -132,7 +132,7 @@ class Energy_Tests(unittest.TestCase):
             xi=xi,
             D=D,
             Instrument=R,
-            Projection=P,
+            Distribution=Dist,
             nonlinearity=f,
             ht=ht,
             N=N,
@@ -160,13 +160,13 @@ class Curvature_Tests(unittest.TestCase):
         ht = ift.HarmonicTransformOperator(hspace, space)
         binbounds = ift.PowerSpace.useful_binbounds(hspace, logarithmic=True)
         pspace = ift.PowerSpace(hspace, binbounds=binbounds)
-        P = ift.PowerProjectionOperator(domain=hspace, power_space=pspace)
+        Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
         xi = ift.Field.from_random(domain=hspace, random_type='normal')
 
         def pspec(k): return 64 / (1 + k**2)**dim
         pspec = ift.PS_field(pspace, pspec)
         tau0 = ift.log(pspec)
-        A = P.adjoint_times(ift.sqrt(pspec))
+        A = Dist(ift.sqrt(pspec))
         n = ift.Field.from_random(domain=space, random_type='normal', std=.01)
         N = ift.DiagonalOperator(n**2)
         s = xi * A
diff --git a/test/test_operators/test_adjoint.py b/test/test_operators/test_adjoint.py
index 0579f857d261c3f7bf4cf63cf318137f92dec043..078b6663c4599f5aee24ffb59ad9bd0871802587 100644
--- a/test/test_operators/test_adjoint.py
+++ b/test/test_operators/test_adjoint.py
@@ -35,15 +35,15 @@ _p_spaces = _p_RG_spaces + [ift.HPSpace(17), ift.GLSpace(8, 13)]
 class Consistency_Tests(unittest.TestCase):
     @expand(product(_h_spaces, [np.float64, np.complex128]))
     def testPPO(self, sp, dtype):
-        op = ift.PowerProjectionOperator(sp)
+        op = ift.PowerDistributor(target=sp)
         ift.extra.consistency_check(op, dtype, dtype)
         ps = ift.PowerSpace(
             sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=False, nbin=3))
-        op = ift.PowerProjectionOperator(sp, ps)
+        op = ift.PowerDistributor(target=sp, power_space=ps)
         ift.extra.consistency_check(op, dtype, dtype)
         ps = ift.PowerSpace(
             sp, ift.PowerSpace.useful_binbounds(sp, logarithmic=True, nbin=3))
-        op = ift.PowerProjectionOperator(sp, ps)
+        op = ift.PowerDistributor(target=sp, power_space=ps)
         ift.extra.consistency_check(op, dtype, dtype)
 
     @expand(product(_h_RG_spaces+_p_RG_spaces,