diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py
index 865938e798b16ca2fa40bb5cd569ace4099cc685..b6bd4f6046fd4ff755523892a7d7f2bdc8e93cff 100644
--- a/demos/log_normal_wiener_filter.py
+++ b/demos/log_normal_wiener_filter.py
@@ -66,19 +66,18 @@ if __name__ == "__main__":
     # m3 = fft(me3[0].position)
 
     # Plotting
-    ift.plotting.plot(mock_signal.real, name='mock_signal.png',
-                      colormap="plasma", xlabel="Pixel Index",
-                      ylabel="Pixel Index")
-    logdata = np.log(ift.dobj.to_global_data(data.val.real)).reshape(signal_space.shape)
+    ift.plotting.plot(mock_signal, name='mock_signal.png', colormap="plasma",
+                      xlabel="Pixel Index", ylabel="Pixel Index")
+    logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
     ift.plotting.plot(ift.Field(signal_space,
                                 val=ift.dobj.from_global_data(logdata)),
                       name="log_of_data.png", colormap="plasma",
                       xlabel="Pixel Index", ylabel="Pixel Index")
-    # ift.plotting.plot(m1.real,name='m_LBFGS.png', colormap="plasma",
+    # ift.plotting.plot(m1,name='m_LBFGS.png', colormap="plasma",
     #                   xlabel="Pixel Index", ylabel="Pixel Index")
-    ift.plotting.plot(m2.real, name='m_Newton.png', colormap="plasma",
+    ift.plotting.plot(m2, name='m_Newton.png', colormap="plasma",
                       xlabel="Pixel Index", ylabel="Pixel Index")
-    # ift.plotting.plot(m3.real, name='m_SteepestDescent.png',
+    # ift.plotting.plot(m3, name='m_SteepestDescent.png',
     #                   colormap="plasma", xlabel="Pixel Index",
     #                   ylabel="Pixel Index")
 
diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py
index d547c52d97e2d2bc8086f66dcda6538350a5a8ee..ae0059575a2a94a8b814f61ea499d439dd81a84b 100644
--- a/demos/nonlinear_critical_filter.py
+++ b/demos/nonlinear_critical_filter.py
@@ -46,7 +46,7 @@ if __name__ == "__main__":
     # Instrument = SmoothingOperator(s_space, sigma=0.01)
     mask = np.ones(s_space.shape)
     mask[6000:8000] = 0.
-    mask = ift.Field(s_space, val=mask)
+    mask = ift.Field(s_space, val=ift.dobj.from_global_data(mask))
 
     MaskOperator = ift.DiagonalOperator(mask)
     InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py
index e32f4836223d8745221bf3473fc0f21208d11cb8..030ffb320d3a9703f7071ce96179c5c0166ed261 100644
--- a/demos/paper_demos/cartesian_wiener_filter.py
+++ b/demos/paper_demos/cartesian_wiener_filter.py
@@ -53,7 +53,6 @@ if __name__ == "__main__":
                      ift.dobj.to_global_data(mock_power_2.val))))
 
     diagonal = ift.power_synthesize_special(mock_power, spaces=(0, 1))**2
-    diagonal = diagonal.real
 
     S = ift.DiagonalOperator(diagonal)
 
diff --git a/demos/wiener_filter_easy.py b/demos/wiener_filter_easy.py
index 2ff3709e729c87a7b2867f4c020a016999743faa..7fe7fb08c4780a7aa286a58323a0be4cb9bad26b 100644
--- a/demos/wiener_filter_easy.py
+++ b/demos/wiener_filter_easy.py
@@ -64,7 +64,7 @@ if __name__ == "__main__":
 
     Sh = ift.create_power_operator(h_space, power_spectrum=pow_spec)
 
-    sp = ift.Field(p_space, val=pow_spec(p_space.k_lengths))
+    sp = ift.PS_field(p_space, pow_spec)
     sh = ift.power_synthesize(sp, real_signal=True)
     ss = fft.inverse_times(sh)
 
diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py
index 3758085d47a13b05720d79ad5a263c75cb367d8e..d79c837f910a1b375473613b34849b2af1c46958 100644
--- a/demos/wiener_filter_via_curvature.py
+++ b/demos/wiener_filter_via_curvature.py
@@ -47,7 +47,6 @@ if __name__ == "__main__":
 
     mock_power = ift.PS_field(power_space, power_spectrum)
     mock_harmonic = ift.power_synthesize(mock_power, real_signal=True)
-    mock_harmonic = mock_harmonic.real
     mock_signal = fft(mock_harmonic)
 
     exposure = 1.
@@ -78,9 +77,9 @@ if __name__ == "__main__":
 
     sspace2 = ift.RGSpace(shape, distances=L/N_pixels/nu.m)
 
-    ift.plotting.plot(ift.Field(sspace2, mock_signal.real.val)/nu.K,
+    ift.plotting.plot(ift.Field(sspace2, mock_signal.val)/nu.K,
                       name="mock_signal.png")
-    data = ift.dobj.to_global_data(data.val.real).reshape(sspace2.shape)/nu.K
+    data = ift.dobj.to_global_data(data.val).reshape(sspace2.shape)/nu.K
     data = ift.Field(sspace2, val=ift.dobj.from_global_data(data))/nu.K
     ift.plotting.plot(ift.Field(sspace2, val=data), name="data.png")
-    ift.plotting.plot(ift.Field(sspace2, m_s.real.val)/nu.K, name="map.png")
+    ift.plotting.plot(ift.Field(sspace2, m_s.val)/nu.K, name="map.png")
diff --git a/nifty/data_objects/random.py b/nifty/data_objects/random.py
index 15d3cbda9100107c450c63ba42b50dea976fdd91..76c222fe8b35a0df4022ad56aa217c2c9c6dfebe 100644
--- a/nifty/data_objects/random.py
+++ b/nifty/data_objects/random.py
@@ -23,7 +23,7 @@ import numpy as np
 class Random(object):
     @staticmethod
     def pm1(dtype, shape):
-        if issubclass(dtype, (complex, np.complexfloating)):
+        if np.issubdtype(dtype, np.complexfloating):
             x = np.array([1+0j, 0+1j, -1+0j, 0-1j], dtype=dtype)
             x = x[np.random.randint(4, size=shape)]
         else:
@@ -32,7 +32,7 @@ class Random(object):
 
     @staticmethod
     def normal(dtype, shape, mean=0., std=1.):
-        if issubclass(dtype, (complex, np.complexfloating)):
+        if np.issubdtype(dtype, np.complexfloating):
             x = np.empty(shape, dtype=dtype)
             x.real = np.random.normal(mean.real, std*np.sqrt(0.5), shape)
             x.imag = np.random.normal(mean.imag, std*np.sqrt(0.5), shape)
@@ -42,12 +42,11 @@ class Random(object):
 
     @staticmethod
     def uniform(dtype, shape, low=0., high=1.):
-        if issubclass(dtype, (complex, np.complexfloating)):
+        if np.issubdtype(dtype, np.complexfloating):
             x = np.empty(shape, dtype=dtype)
             x.real = np.random.uniform(low, high, shape)
             x.imag = np.random.uniform(low, high, shape)
-        elif dtype in [np.dtype('int8'), np.dtype('int16'), np.dtype('int32'),
-                       np.dtype('int64')]:
+        elif np.issubdtype(dtype, np.integer):
             x = np.random.random.randint(low, high+1, shape)
         else:
             x = np.random.uniform(low, high, shape)
diff --git a/nifty/field.py b/nifty/field.py
index 03fe28e9e5e67d947b8204ab6e7028e63f9d466d..c80d75a49f3c2f86b46eff54a99d6be80f18d386 100644
--- a/nifty/field.py
+++ b/nifty/field.py
@@ -213,11 +213,15 @@ class Field(object):
     @property
     def real(self):
         """ The real part of the field (data is not copied)."""
+        if not np.issubdtype(self.dtype, np.complexfloating):
+            raise ValueError(".real called on a non-complex Field")
         return Field(self.domain, self.val.real)
 
     @property
     def imag(self):
         """ The imaginary part of the field (data is not copied)."""
+        if not np.issubdtype(self.dtype, np.complexfloating):
+            raise ValueError(".imag called on a non-complex Field")
         return Field(self.domain, self.val.imag)
 
     def copy(self):
diff --git a/nifty/library/critical_power_energy.py b/nifty/library/critical_power_energy.py
index d38739780854ddde7426d6752eb8e0a0c4d321d3..6d007fb1079f06037e01b23999dc6138297b141f 100644
--- a/nifty/library/critical_power_energy.py
+++ b/nifty/library/critical_power_energy.py
@@ -89,7 +89,7 @@ class CriticalPowerEnergy(Energy):
         gradient = -self._theta
         gradient += self.alpha-0.5
         gradient += Tt
-        self._gradient = gradient.real
+        self._gradient = gradient
 
     def at(self, position):
         return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
diff --git a/nifty/operators/diagonal_operator.py b/nifty/operators/diagonal_operator.py
index 44440114df3d1c723cc24ef05efb1ad0756e5c8e..51a9d456b2932667ddb16947a3af0b103020a135 100644
--- a/nifty/operators/diagonal_operator.py
+++ b/nifty/operators/diagonal_operator.py
@@ -140,7 +140,7 @@ class DiagonalOperator(EndomorphicOperator):
     @property
     def self_adjoint(self):
         if self._self_adjoint is None:
-            if not issubclass(self._diagonal.dtype.type, np.complexfloating):
+            if not np.issubdtype(self._diagonal.dtype, np.complexfloating):
                 self._self_adjoint = True
             else:
                 self._self_adjoint = (self._diagonal.val.imag == 0).all()
diff --git a/nifty/operators/fft_operator.py b/nifty/operators/fft_operator.py
index c34ecaac697bf62c77e407d2518eaeb53ff6a9d8..78d13da6f5e52d004e1a86f66fd79a77c15dd2ce 100644
--- a/nifty/operators/fft_operator.py
+++ b/nifty/operators/fft_operator.py
@@ -107,7 +107,7 @@ class FFTOperator(LinearOperator):
             self._trafo = SphericalTransformation(pdom, hdom, space)
 
     def _times_helper(self, x):
-        if issubclass(x.dtype.type, np.complexfloating):
+        if np.issubdtype(x.dtype, np.complexfloating):
             res = (self._trafo.transform(x.real) +
                    1j * self._trafo.transform(x.imag))
         else:
diff --git a/nifty/plotting/plot.py b/nifty/plotting/plot.py
index ce65d23d2868aed158ddea0ee4a1de0a618ca892..88accb5f84acab7e48db65cd309131960d8fd552 100644
--- a/nifty/plotting/plot.py
+++ b/nifty/plotting/plot.py
@@ -45,9 +45,11 @@ def _find_closest(A, target):
 def _makeplot(name):
     import matplotlib.pyplot as plt
     if dobj.rank != 0:
+        plt.close()
         return
     if name is None:
         plt.show()
+        plt.close()
         return
     extension = os.path.splitext(name)[1]
     if extension == ".pdf":
diff --git a/nifty/sugar.py b/nifty/sugar.py
index da15940964f76a90e57c45bede6b019e3a538ffa..23999493b7d4e664f1762d721b45a918dfe80bd6 100644
--- a/nifty/sugar.py
+++ b/nifty/sugar.py
@@ -94,10 +94,17 @@ def power_analyze(field, spaces=None, binbounds=None,
     if len(spaces) == 0:
         raise ValueError("No space for analysis specified.")
 
+    field_real = not np.issubdtype(field.dtype, np.complexfloating)
+    if (not field_real) and keep_phase_information:
+        raise ValueError("cannot keep phase from real-valued input Field")
+
     if keep_phase_information:
         parts = [field.real*field.real, field.imag*field.imag]
     else:
-        parts = [field.real*field.real + field.imag*field.imag]
+        if field_real:
+            parts = [field**2]
+        else:
+            parts = [field.real*field.real + field.imag*field.imag]
 
     parts = [part.weight(1, spaces) for part in parts]
     for space_index in spaces:
@@ -162,6 +169,9 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
     """
 
     spec = _compute_spec(field, spaces)
+    self_real = not np.issubdtype(spec.dtype, np.complexfloating)
+    if (not real_power) and self_real:
+        raise ValueError("can't draw complex realizations from real spectrum")
 
     # create random samples: one or two, depending on whether the
     # power spectrum is real or complex
@@ -176,7 +186,7 @@ def power_synthesize(field, spaces=None, real_power=True, real_signal=True):
         field.from_random('normal', mean=0., std=1.,
                           domain=spec.domain, dtype=np.float)
 
-    result[0] *= spec.real
+    result[0] *= spec if self_real else spec.real
     if not real_power:
         result[1] *= spec.imag
 
@@ -190,7 +200,7 @@ def power_synthesize_special(field, spaces=None):
     field.from_random('normal', mean=0., std=1.,
                       domain=spec.domain, dtype=np.complex)
 
-    return spec.real
+    return spec
 
 
 def create_power_field(domain, power_spectrum, dtype=None):
@@ -206,12 +216,8 @@ def create_power_field(domain, power_spectrum, dtype=None):
     else:
         power_domain = PowerSpace(domain)
         fp = PS_field(power_domain, power_spectrum, dtype)
-    f = PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
-
-    if not issubclass(fp.dtype.type, np.complexfloating):
-        f = f.real
 
-    return f
+    return PowerProjectionOperator(domain, power_domain).adjoint_times(fp)
 
 
 def create_power_operator(domain, power_spectrum, space=None, dtype=None):
diff --git a/nifty/utilities.py b/nifty/utilities.py
index ac642f33e621f50ad2000feee67339fc4df8c303..107e725fcfed24a6168be1ac6b99689bf8be27f6 100644
--- a/nifty/utilities.py
+++ b/nifty/utilities.py
@@ -129,7 +129,7 @@ def hartley(a, axes=None):
     if axes is not None and \
             not all(axis < len(a.shape) for axis in axes):
         raise ValueError("Provided axes do not match array shape")
-    if issubclass(a.dtype.type, np.complexfloating):
+    if np.issubdtype(a.dtype, np.complexfloating):
         raise TypeError("Hartley transform requires real-valued arrays.")
 
     from pyfftw.interfaces.numpy_fft import rfftn
@@ -171,7 +171,7 @@ def my_fftn_r2c(a, axes=None):
     if axes is not None and \
             not all(axis < len(a.shape) for axis in axes):
         raise ValueError("Provided axes do not match array shape")
-    if issubclass(a.dtype.type, np.complexfloating):
+    if np.issubdtype(a.dtype, np.complexfloating):
         raise TypeError("Transform requires real-valued input arrays.")
 
     from pyfftw.interfaces.numpy_fft import rfftn