diff --git a/gaussian_processes.py b/gaussian_processes.py
new file mode 100644
index 0000000000000000000000000000000000000000..998ee07c3ee722a17571947e52fb5116412174e4
--- /dev/null
+++ b/gaussian_processes.py
@@ -0,0 +1,59 @@
+import nifty8 as ift
+
+position_space = ift.RGSpace((128,128))
+harmonic_space = position_space.get_default_codomain()
+HT = ift.HarmonicTransformOperator(harmonic_space)
+print(HT)
+print(HT.domain, HT.target)
+
+power_space = ift.PowerSpace(harmonic_space)
+power_spectrum = (lambda k: 1./(10.+ k**4))
+power_spectrum = ift.PS_field(power_space, power_spectrum)
+
+pl = ift.Plot()
+pl.add(power_spectrum)
+pl.output()
+
+PD = ift.PowerDistributor(harmonic_space)
+spectrum = PD(power_spectrum)
+
+xi = ift.FieldAdapter(harmonic_space, 'inp')
+
+# s = HT(sqrt(power_spectrum) * xi)
+s = HT @ ift.DiagonalOperator(spectrum.sqrt()) @ xi
+# s = s.exp()
+print(s)
+print(s.domain, s.target)
+
+pl = ift.Plot()
+for _ in range(6):
+    pl.add(s(ift.from_random(s.domain)))
+pl.output()
+
+args = {
+    'offset_mean': 0,
+    'offset_std': (1e-3, 1e-6),
+
+    # Amplitude of field fluctuations
+    'fluctuations': (1., 0.8),  # 1.0, 1e-2
+
+    # Exponent of power law power spectrum component
+    'loglogavgslope': (-3., 1),  # -6.0, 1
+
+    # Amplitude of integrated Wiener process power spectrum component
+    'flexibility': (2, 1.),  # 1.0, 0.5
+
+    # How ragged the integrated Wiener process component is
+    'asperity': (0.2, 0.1)  # 0.1, 0.5
+}
+correlated_field = ift.SimpleCorrelatedField(position_space, **args)
+power_spectrum = correlated_field.power_spectrum
+
+pl = ift.Plot()
+pl.add([power_spectrum(ift.from_random(power_spectrum.domain)) for _ in range(10)])
+pl.output()
+
+pl = ift.Plot()
+for _ in range(6):
+    pl.add(correlated_field(ift.from_random(correlated_field.domain)))
+pl.output()
\ No newline at end of file