diff --git a/docs/Makefile b/docs/Makefile
deleted file mode 100644
index 53bfc172ec19088227436bd75c551eae5d9e321e..0000000000000000000000000000000000000000
--- a/docs/Makefile
+++ /dev/null
@@ -1,216 +0,0 @@
-# Makefile for Sphinx documentation
-#
-
-# You can set these variables from the command line.
-SPHINXOPTS    =
-SPHINXBUILD   = sphinx-build
-PAPER         =
-BUILDDIR      = build
-
-# User-friendly check for sphinx-build
-ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1)
-$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don't have Sphinx installed, grab it from http://sphinx-doc.org/)
-endif
-
-# Internal variables.
-PAPEROPT_a4     = -D latex_paper_size=a4
-PAPEROPT_letter = -D latex_paper_size=letter
-ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
-# the i18n builder cannot share the environment and doctrees with the others
-I18NSPHINXOPTS  = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
-
-.PHONY: help
-help:
-	@echo "Please use \`make <target>' where <target> is one of"
-	@echo "  html       to make standalone HTML files"
-	@echo "  dirhtml    to make HTML files named index.html in directories"
-	@echo "  singlehtml to make a single large HTML file"
-	@echo "  pickle     to make pickle files"
-	@echo "  json       to make JSON files"
-	@echo "  htmlhelp   to make HTML files and a HTML help project"
-	@echo "  qthelp     to make HTML files and a qthelp project"
-	@echo "  applehelp  to make an Apple Help Book"
-	@echo "  devhelp    to make HTML files and a Devhelp project"
-	@echo "  epub       to make an epub"
-	@echo "  latex      to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
-	@echo "  latexpdf   to make LaTeX files and run them through pdflatex"
-	@echo "  latexpdfja to make LaTeX files and run them through platex/dvipdfmx"
-	@echo "  text       to make text files"
-	@echo "  man        to make manual pages"
-	@echo "  texinfo    to make Texinfo files"
-	@echo "  info       to make Texinfo files and run them through makeinfo"
-	@echo "  gettext    to make PO message catalogs"
-	@echo "  changes    to make an overview of all changed/added/deprecated items"
-	@echo "  xml        to make Docutils-native XML files"
-	@echo "  pseudoxml  to make pseudoxml-XML files for display purposes"
-	@echo "  linkcheck  to check all external links for integrity"
-	@echo "  doctest    to run all doctests embedded in the documentation (if enabled)"
-	@echo "  coverage   to run coverage check of the documentation (if enabled)"
-
-.PHONY: clean
-clean:
-	rm -rf $(BUILDDIR)/*
-
-.PHONY: html
-html:
-	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
-	@echo
-	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
-
-.PHONY: dirhtml
-dirhtml:
-	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
-	@echo
-	@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
-
-.PHONY: singlehtml
-singlehtml:
-	$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
-	@echo
-	@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
-
-.PHONY: pickle
-pickle:
-	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
-	@echo
-	@echo "Build finished; now you can process the pickle files."
-
-.PHONY: json
-json:
-	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
-	@echo
-	@echo "Build finished; now you can process the JSON files."
-
-.PHONY: htmlhelp
-htmlhelp:
-	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
-	@echo
-	@echo "Build finished; now you can run HTML Help Workshop with the" \
-	      ".hhp project file in $(BUILDDIR)/htmlhelp."
-
-.PHONY: qthelp
-qthelp:
-	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
-	@echo
-	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
-	      ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
-	@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/NIFTY.qhcp"
-	@echo "To view the help file:"
-	@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/NIFTY.qhc"
-
-.PHONY: applehelp
-applehelp:
-	$(SPHINXBUILD) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp
-	@echo
-	@echo "Build finished. The help book is in $(BUILDDIR)/applehelp."
-	@echo "N.B. You won't be able to view it unless you put it in" \
-	      "~/Library/Documentation/Help or install it in your application" \
-	      "bundle."
-
-.PHONY: devhelp
-devhelp:
-	$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
-	@echo
-	@echo "Build finished."
-	@echo "To view the help file:"
-	@echo "# mkdir -p $$HOME/.local/share/devhelp/NIFTY"
-	@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/NIFTY"
-	@echo "# devhelp"
-
-.PHONY: epub
-epub:
-	$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
-	@echo
-	@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
-
-.PHONY: latex
-latex:
-	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
-	@echo
-	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
-	@echo "Run \`make' in that directory to run these through (pdf)latex" \
-	      "(use \`make latexpdf' here to do that automatically)."
-
-.PHONY: latexpdf
-latexpdf:
-	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
-	@echo "Running LaTeX files through pdflatex..."
-	$(MAKE) -C $(BUILDDIR)/latex all-pdf
-	@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
-
-.PHONY: latexpdfja
-latexpdfja:
-	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
-	@echo "Running LaTeX files through platex and dvipdfmx..."
-	$(MAKE) -C $(BUILDDIR)/latex all-pdf-ja
-	@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
-
-.PHONY: text
-text:
-	$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
-	@echo
-	@echo "Build finished. The text files are in $(BUILDDIR)/text."
-
-.PHONY: man
-man:
-	$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
-	@echo
-	@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
-
-.PHONY: texinfo
-texinfo:
-	$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
-	@echo
-	@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
-	@echo "Run \`make' in that directory to run these through makeinfo" \
-	      "(use \`make info' here to do that automatically)."
-
-.PHONY: info
-info:
-	$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
-	@echo "Running Texinfo files through makeinfo..."
-	make -C $(BUILDDIR)/texinfo info
-	@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
-
-.PHONY: gettext
-gettext:
-	$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
-	@echo
-	@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
-
-.PHONY: changes
-changes:
-	$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
-	@echo
-	@echo "The overview file is in $(BUILDDIR)/changes."
-
-.PHONY: linkcheck
-linkcheck:
-	$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
-	@echo
-	@echo "Link check complete; look for any errors in the above output " \
-	      "or in $(BUILDDIR)/linkcheck/output.txt."
-
-.PHONY: doctest
-doctest:
-	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
-	@echo "Testing of doctests in the sources finished, look at the " \
-	      "results in $(BUILDDIR)/doctest/output.txt."
-
-.PHONY: coverage
-coverage:
-	$(SPHINXBUILD) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage
-	@echo "Testing of coverage in the sources finished, look at the " \
-	      "results in $(BUILDDIR)/coverage/python.txt."
-
-.PHONY: xml
-xml:
-	$(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml
-	@echo
-	@echo "Build finished. The XML files are in $(BUILDDIR)/xml."
-
-.PHONY: pseudoxml
-pseudoxml:
-	$(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml
-	@echo
-	@echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml."
diff --git a/docs/howtogenerate.txt b/docs/howtogenerate.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c0a2cdb60377735229c039f05615268c67b71af6
--- /dev/null
+++ b/docs/howtogenerate.txt
@@ -0,0 +1,10 @@
+Two step creation of webpages:
+
+sphinx-apidoc -l -e -d 3 -o sphinx/source/mod/ nifty/ nifty/plotting/ nifty/spaces/power_space/power_indices.py nifty/spaces/power_space/power_index_factory.py nifty/config/ nifty/basic_arithmetics.py nifty/nifty_meta.py nifty/random.py nifty/version.py  nifty/field_types/ nifty/operators/fft_operator/transformations/rg_transforms.py
+
+
+creates all .rst files neccesary for ModuleIndex excluding helper modules
+
+sphinx-build -b html sphinx/source/ sphinx/build/
+
+generates html filel amd build directory
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 8b422952541ea795c53fe8a3c62ddb8098f811ae..a68fcf8d3bacaebb26d759b94dc5cb87a09b6c79 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -16,6 +16,9 @@ from nifty import *
 import sys
 import os
 
+import sphinx_rtd_theme
+
+
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
@@ -77,7 +80,7 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'NIFTY'
-copyright = u'2017, Theo Steininger'
+copyright = u'2013-2017, Max-Planck-Society'
 author = u'Theo Steininger'
 
 # The version info for the project you're documenting, acts as replacement for
@@ -85,9 +88,9 @@ author = u'Theo Steininger'
 # built documents.
 #
 # The short X.Y version.
-version = u'3.0.4'
+version = u'3.0'
 # The full version, including alpha/beta/rc tags.
-release = u'3.0.x'
+release = u'3.0.4'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -138,18 +141,19 @@ todo_include_todos = True
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-html_theme = 'alabaster'
-#html_theme = 'classic'
+html_theme = "sphinx_rtd_theme"
 
+# Theme options are theme-specific and customize the look and feel of a theme
+# further.  For a list of options available for each theme, see the
+# documentation.
 html_theme_options = {
-    'page_width':''
+     'collapse_navigation': False,
+     'display_version': False,
+     'navigation_depth': 3,
 }
 
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
-# Theme options are theme-specific and customize the look and feel of a theme
-# further.  For a list of options available for each theme, see the
-# documentation.
-#html_theme_options = {}
 
 # Add any paths that contain custom themes here, relative to this directory.
 #html_theme_path = []
@@ -163,7 +167,7 @@ html_theme_options = {
 
 # The name of an image file (relative to this directory) to place at the top
 # of the sidebar.
-html_logo = 'nifty_logo_black.png'
+html_logo = 'nifty_logo_white.png'
 
 # The name of an image file (relative to this directory) to use as a favicon of
 # the docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
@@ -196,7 +200,7 @@ html_last_updated_fmt = '%b %d, %Y'
 #html_additional_pages = {}
 
 # If false, no module index is generated.
-#html_domain_indices = True
+html_domain_indices = False
 
 # If false, no index is generated.
 #html_use_index = True
diff --git a/docs/source/diagonal_operator.rst b/docs/source/diagonal_operator.rst
deleted file mode 100644
index 0dc633e47396c4349809c813f9c7d7cbe6e737ab..0000000000000000000000000000000000000000
--- a/docs/source/diagonal_operator.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``DiagonalOperator`` class -- ...
-.....................................
-
-.. autoclass:: DiagonalOperator
-    :show-inheritance:
-    :members:
diff --git a/docs/source/fft_operator.rst b/docs/source/fft_operator.rst
deleted file mode 100644
index 9bc6dc205473981468cfb27a46d0b8d3381c54ef..0000000000000000000000000000000000000000
--- a/docs/source/fft_operator.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``FFTOperator`` class -- Fourier Transformations
-....................................................
-
-.. autoclass:: FFTOperator
-    :show-inheritance:
-    :members:
diff --git a/docs/source/field.rst b/docs/source/field.rst
index fe75f39ff67c62161de1babe77fdc7b6826a6543..eb70e5e514586f1bc7274b54112bb24f08dd8d95 100644
--- a/docs/source/field.rst
+++ b/docs/source/field.rst
@@ -11,7 +11,3 @@ In NIFTY, Fields are used to store data arrays and carry all the needed metainfo
     :show-inheritance:
     :members:
 
-    .. rubric:: Methods
-
-    .. autoautosummary:: Field
-        :methods:
diff --git a/docs/source/gallery.rst b/docs/source/gallery.rst
new file mode 100644
index 0000000000000000000000000000000000000000..2d556fc21c4b8318e8b5da696dff27e759c544ef
--- /dev/null
+++ b/docs/source/gallery.rst
@@ -0,0 +1,137 @@
+Image Gallery
+-------------
+
+Transformations & Projections
+.............................
+
+.. currentmodule:: nifty
+
+The "Faraday Map" [1]_ in spherical representation on a :py:class:`hp_space` and a :py:class:`gl_space`, their quadrupole projections, the uncertainty of the map, and the angular power spectrum.
+
++----------------------------+----------------------------+
+| .. image:: images/f_00.png | .. image:: images/f_01.png |
+|     :width:  90 %          |     :width:  90 %          |
++----------------------------+----------------------------+
+| .. image:: images/f_02.png | .. image:: images/f_03.png |
+|     :width:  90 %          |     :width:  90 %          |
++----------------------------+----------------------------+
+| .. image:: images/f_04.png | .. image:: images/f_05.png |
+|     :width:  90 %          |     :width:  70 %          |
++----------------------------+----------------------------+
+
+Gaussian random fields
+......................
+
+Statistically homogeneous and isotropic Gaussian random fields drawn from different power spectra.
+
++----------------------------+----------------------------+
+| .. image:: images/t_03.png | .. image:: images/t_04.png |
+|     :width:  60 %          |     :width:  70 %          |
++----------------------------+----------------------------+
+| .. image:: images/t_05.png | .. image:: images/t_06.png |
+|     :width:  60 %          |     :width:  70 %          |
++----------------------------+----------------------------+
+
+
+Wiener filtering I
+..................
+
+Wiener filter reconstruction of Gaussian random signal.
+
++--------------------------------+--------------------------------+--------------------------------+
+| original signal                | noisy data                     | reconstruction                 |
++================================+================================+================================+
+| .. image:: images/rg1_s.png    | .. image:: images/rg1_d.png    | .. image:: images/rg1_m.png    |
+|     :width:  90 %              |     :width:  90 %              |     :width:  90 %              |
++--------------------------------+--------------------------------+--------------------------------+
+| .. image:: images/rg2_s_pm.png | .. image:: images/rg2_d_pm.png | .. image:: images/rg2_m_pm.png |
+|     :width:  90 %              |     :width:  90 %              |     :width:  90 %              |
++--------------------------------+--------------------------------+--------------------------------+
+| .. image:: images/hp_s.png     | .. image:: images/hp_d.png     | .. image:: images/hp_m.png     |
+|     :width:  90 %              |     :width:  90 %              |     :width:  90 %              |
++--------------------------------+--------------------------------+--------------------------------+
+
+Image reconstruction
+....................
+
+Image reconstruction of the classic "Moon Surface" image. The original image "Moon Surface" was taken from the `USC-SIPI image database <http://sipi.usc.edu/database/>`_.
+
++-----------------------------------+-----------------------------------+-----------------------------------+
+| .. image:: images/moon_s.png      | .. image:: images/moon_d.png      | .. image:: images/moon_m.png      |
+|     :width:  90 %                 |     :width:  90 %                 |     :width:  90 %                 |
++-----------------------------------+-----------------------------------+-----------------------------------+
+| .. image:: images/moon_kernel.png | .. image:: images/moon_mask.png   | .. image:: images/moon_sigma.png  |
+|     :width:  90 %                 |     :width:  90 %                 |     :width:  90 %                 |
++-----------------------------------+-----------------------------------+-----------------------------------+
+
+Wiener filtering II
+...................
+
+Wiener filter reconstruction results for the full and partially blinded data. Shown are the original signal (orange), the reconstruction (green), and :math:`1\sigma`-confidence interval (gray).
+
++--------------------------------------+--------------------------------------+
+| noisy data                           | reconstruction results               |
++======================================+======================================+
+| .. image:: images/rg1_d.png          | .. image:: images/rg1_m_err_.png     |
+|     :width:  90 %                    |     :width:  90 %                    |
++--------------------------------------+--------------------------------------+
+| .. image:: images/rg1_d_gap.png      | .. image:: images/rg1_m_gap_err_.png |
+|     :width:  90 %                    |     :width:  90 %                    |
++--------------------------------------+--------------------------------------+
+
+D\ :sup:`3`\ PO -- Denoising, Deconvolving, and Decomposing Photon Observations
+...............................................................................
+
+Application of the D\ :sup:`3`\ PO algorithm [2]_ showing the raw photon count data and the denoised, deconvolved, and decomposed reconstruction of the diffuse photon flux.
+
++--------------------------------------+--------------------------------------+
+| .. image:: images/D3PO_data.png      | .. image:: images/D3PO_diffuse.png   |
+|     :width:  95 %                    |     :width:  95 %                    |
++--------------------------------------+--------------------------------------+
+
+RESOLVE -- Aperature synthesis imaging in radio astronomy
+.........................................................
+
+Signal inference on simulated single-frequency data: reconstruction by CLEAN (using uniform weighting) and by RESOLVE [3]_ (using IFT & NIFTY).
+
++-------------------------------------+-------------------------------------+-------------------------------------+
+| .. image:: images/radio_signal.png  | .. image:: images/radio_CLEAN.png   | .. image:: images/radio_RESOLVE.png |
+|     :width:  90 %                   |     :width:  90 %                   |     :width:  90 %                   |
++-------------------------------------+-------------------------------------+-------------------------------------+
+
+D\ :sup:`3`\ PO -- light
+........................
+
+Inference of the mock distribution of some species across Australia exploiting geospatial correlations in a (strongly) simplified scenario [4]_.
+
++--------------------------------+--------------------------------+--------------------------------+
+| .. image:: images/au_data.png  | .. image:: images/au_map.png   | .. image:: images/au_error.png |
+|     :width:  90 %              |     :width:  90 %              |     :width:  90 %              |
++--------------------------------+--------------------------------+--------------------------------+
+
+NIFTY meets Lensing
+...................
+
+Signal reconstruction for a simulated image that has undergone strong gravitational lensing. Without *a priori* knowledge of the signal covariance :math:`S`, a common approach rescaling the  `Laplace-Operator <http://de.wikipedia.org/wiki/Laplace-Operator>`_ and IFT's `"critical" filter <./demo_excaliwir.html#critical-wiener-filtering>`_ are compared.
+
++--------------------------------+--------------------------------+--------------------------------+--------------------------------+
+| .. image:: images/lens_s0.png  | .. image:: images/lens_d0.png  | .. image:: images/lens_m1.png  | .. image:: images/lens_m2.png  |
+|     :width:  80 %              |     :width:  80 %              |     :width:  80 %              |     :width:  80 %              |
+|                                |                                |                                |                                |
+|                                |                                | .. math::                      | .. math::                      |
+|                                |                                |     S(x,y) &=                  |     S(x,y) &=                  |
+|                                |                                |         \lambda \: \Delta^{-1} |         S(|x-y|)               |
+|                                |                                |         \\ \equiv              |         \\ \equiv              |
+|                                |                                |     S(k,l) &= \delta(k-l)      |     S(k,l) &= \delta(k-l)      |
+|                                |                                |         \: \lambda \: k^{-2}   |         \: P(k)                |
++--------------------------------+--------------------------------+--------------------------------+--------------------------------+
+
+.. [1] N. Oppermann et. al., "An improved map of the Galactic Faraday sky", Astronomy & Astrophysics, vol. 542, id. A93, p. 14, see also the `project homepage <http://www.mpa-garching.mpg.de/ift/faraday/>`_
+
+.. [2] M. Selig et. al., "Denoising, Deconvolving, and Decomposing Photon Observations", submitted to Astronomy & Astrophysics, 2013; `arXiv:1311.1888 <http://www.arxiv.org/abs/1311.1888>`_
+
+.. [3] H. Junklewitz et. al., "RESOLVE: A new algorithm for aperture synthesis imaging of extended emission in radio astronomy", submitted to Astronomy & Astrophysics, 2013; `arXiv:1311.5282 <http://www.arxiv.org/abs/1311.5282>`_
+
+.. [4] M. Selig, "The NIFTY way of Bayesian signal inference", submitted proceeding of the 33rd MaxEnt, 2013
+
+
diff --git a/docs/source/ift.rst b/docs/source/ift.rst
new file mode 100644
index 0000000000000000000000000000000000000000..0f61d5cbae7266401843b67162d47925a1e22676
--- /dev/null
+++ b/docs/source/ift.rst
@@ -0,0 +1,83 @@
+IFT -- Information Field Theory
+===============================
+
+Theoretical Background
+----------------------
+
+
+`Information Field Theory <http://www.mpa-garching.mpg.de/ift/>`_ [1]_ (IFT) is information theory, the logic of reasoning under uncertainty, applied to fields. A field can be any quantity defined over some space, e.g. the air temperature over Europe, the magnetic field strength in the Milky Way, or the matter density in the Universe. IFT describes how data and knowledge can be used to infer field properties. Mathematically it is a statistical field theory and exploits many of the tools developed for such. Practically, it is a framework for signal processing and image reconstruction.
+
+IFT is fully Bayesian. How else can infinitely many field degrees of freedom be constrained by finite data?
+
+It can be used without the knowledge of Feynman diagrams. There is a full toolbox of methods. It reproduces many known well working algorithms. This should be reassuring. And, there were certainly previous works in a similar spirit. Anyhow, in many cases IFT provides novel rigorous ways to extract information from data.
+
+.. tip:: An *in-a-nutshell introduction to information field theory* can be found in [2]_.
+
+.. [1] T. Ensslin et al., "Information field theory for cosmological perturbation reconstruction and nonlinear signal analysis", PhysRevD.80.105005, 09/2009; `arXiv:0806.3474 <http://www.arxiv.org/abs/0806.3474>`_
+
+.. [2] T. Ensslin, "Information field theory", accepted for the proceedings of MaxEnt 2012 -- the 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering; `arXiv:1301.2556 <http://arxiv.org/abs/1301.2556>`_
+
+
+Discretized continuum
+---------------------
+
+The representation of fields that are mathematically defined on a continuous space in a finite computer environment is a common necessity. The goal hereby is to preserve the continuum limit in the calculus in order to ensure a resolution independent discretization.
+
++-----------------------------+-----------------------------+
+| .. image:: images/42vs6.png | .. image:: images/42vs9.png |
+|     :width:  100 %          |     :width:  100 %          |
++-----------------------------+-----------------------------+
+
+Any partition of the continuous position space :math:`\Omega` (with volume :math:`V`) into a set of :math:`Q` disjoint, proper subsets :math:`\Omega_q` (with volumes :math:`V_q`) defines a pixelization,
+
+.. math::
+
+    \Omega &\quad=\quad \dot{\bigcup_q} \; \Omega_q \qquad \mathrm{with} \qquad q \in \{1,\dots,Q\} \subset \mathbb{N}
+    , \\
+    V &\quad=\quad \int_\Omega \mathrm{d}x \quad=\quad \sum_{q=1}^Q \int_{\Omega_q} \mathrm{d}x \quad=\quad \sum_{q=1}^Q V_q
+    .
+
+Here the number :math:`Q` characterizes the resolution of the pixelization and the continuum limit is described by :math:`Q \rightarrow \infty` and :math:`V_q \rightarrow 0` for all :math:`q \in \{1,\dots,Q\}` simultaneously. Moreover, the above equation defines a discretization of continuous integrals, :math:`\int_\Omega \mathrm{d}x \mapsto \sum_q V_q`.
+
+Any valid discretization scheme for a field :math:`{s}` can be described by a mapping,
+
+.. math::
+
+    s(x \in \Omega_q) \quad\mapsto\quad s_q \quad=\quad \int_{\Omega_q} \mathrm{d}x \; w_q(x) \; s(x)
+    ,
+
+if the weighting function :math:`w_q(x)` is chosen appropriately. In order for the discretized version of the field to converge to the actual field in the continuum limit, the weighting functions need to be normalized in each subset; i.e., :math:`\forall q: \int_{\Omega_q} \mathrm{d}x \; w_q(x) = 1`. Choosing such a weighting function that is constant with respect to :math:`x` yields
+
+.. math::
+
+    s_q = \frac{\int_{\Omega_q} \mathrm{d}x \; s(x)}{\int_{\Omega_q} \mathrm{d}x} = \left< s(x) \right>_{\Omega_q}
+    ,
+
+which corresponds to a discretization of the field by spatial averaging. Another common and equally valid choice is :math:`w_q(x) = \delta(x-x_q)`, which distinguishes some position :math:`x_q \in \Omega_q`, and evaluates the continuous field at this position,
+
+.. math::
+
+    s_q \quad=\quad \int_{\Omega_q} \mathrm{d}x \; \delta(x-x_q) \; s(x) \quad=\quad s(x_q)
+    .
+
+In practice, one often makes use of the spatially averaged pixel position, :math:`x_q = \left< x \right>_{\Omega_q}`. If the resolution is high enough to resolve all features of the signal field :math:`{s}`, both of these discretization schemes approximate each other, :math:`\left< s(x) \right>_{\Omega_q} \approx s(\left< x \right>_{\Omega_q})`, since they approximate the continuum limit by construction. (The approximation of :math:`\left< s(x) \right>_{\Omega_q} \approx s(x_q \in \Omega_q)` marks a resolution threshold beyond which further refinement of the discretization reveals no new features; i.e., no new information content of the field :math:`{s}`.)
+
+All operations involving position integrals can be normalized in accordance with the above definitions. For example, the scalar product between two fields :math:`{s}` and :math:`{u}` is defined as
+
+.. math::
+
+    {s}^\dagger {u} \quad=\quad \int_\Omega \mathrm{d}x \; s^*(x) \; u(x) \quad\approx\quad \sum_{q=1}^Q V_q^{\phantom{*}} \; s_q^* \; u_q^{\phantom{*}}
+    ,
+
+where :math:`\dagger` denotes adjunction and :math:`*` complex conjugation. Since the above approximation becomes an equality in the continuum limit, the scalar product is independent of the pixelization scheme and resolution, if the latter is sufficiently high.
+
+The above line of argumentation analogously applies to the discretization of operators. For a linear operator :math:`{A}` acting on some field :math:`{s}` as :math:`{A} {s} = \int_\Omega \mathrm{d}y \; A(x,y) \; s(y)`, a matrix representation discretized with constant weighting functions is given by
+
+.. math::
+
+    A(x \in \Omega_p, y \in \Omega_q) \quad\mapsto\quad A_{pq} \quad=\quad \frac{\iint_{\Omega_p \Omega_q} \mathrm{d}x \, \mathrm{d}y \; A(x,y)}{\iint_{\Omega_p \Omega_q} \mathrm{d}x \, \mathrm{d}y} \quad=\quad \big< \big< A(x,y) \big>_{\Omega_p} \big>_{\Omega_q}
+    .
+
+The proper discretization of spaces, fields, and operators, as well as the normalization of position integrals, is essential for the conservation of the continuum limit. Their consistent implementation in NIFTY allows a pixelization independent coding of algorithms.
+
+
diff --git a/docs/source/images/42vs6.png b/docs/source/images/42vs6.png
new file mode 100644
index 0000000000000000000000000000000000000000..b867e21a0f4b764e577de9d7d894a35be4b60f94
Binary files /dev/null and b/docs/source/images/42vs6.png differ
diff --git a/docs/source/images/42vs9.png b/docs/source/images/42vs9.png
new file mode 100644
index 0000000000000000000000000000000000000000..c49cbef790b2a90e84df007b87ef15ccb4f56804
Binary files /dev/null and b/docs/source/images/42vs9.png differ
diff --git a/docs/source/images/D3PO_data.png b/docs/source/images/D3PO_data.png
new file mode 100644
index 0000000000000000000000000000000000000000..15fdc42d3633514ca06bcf2d7d812bb7d5c2b326
Binary files /dev/null and b/docs/source/images/D3PO_data.png differ
diff --git a/docs/source/images/D3PO_diffuse.png b/docs/source/images/D3PO_diffuse.png
new file mode 100644
index 0000000000000000000000000000000000000000..b96e4a681d95cc5446ff2a9face14f19736621af
Binary files /dev/null and b/docs/source/images/D3PO_diffuse.png differ
diff --git a/docs/source/images/Skl.png b/docs/source/images/Skl.png
new file mode 100644
index 0000000000000000000000000000000000000000..86f967ab770c6b7780c09d87cd439d54113fce69
Binary files /dev/null and b/docs/source/images/Skl.png differ
diff --git a/docs/source/images/Sxy.png b/docs/source/images/Sxy.png
new file mode 100644
index 0000000000000000000000000000000000000000..1118ce81ddf302aff7c2ed0141c87f2a0d306f6a
Binary files /dev/null and b/docs/source/images/Sxy.png differ
diff --git a/docs/source/images/au_data.png b/docs/source/images/au_data.png
new file mode 100644
index 0000000000000000000000000000000000000000..fe1b778f5c395845c0de3ac8622b1db98fdb51d0
Binary files /dev/null and b/docs/source/images/au_data.png differ
diff --git a/docs/source/images/au_error.png b/docs/source/images/au_error.png
new file mode 100644
index 0000000000000000000000000000000000000000..183bd93c93abe84396d706b4d7dde6b7c2db6659
Binary files /dev/null and b/docs/source/images/au_error.png differ
diff --git a/docs/source/images/au_map.png b/docs/source/images/au_map.png
new file mode 100644
index 0000000000000000000000000000000000000000..32ec5f5884736bb306dfe285c77eeac4da9caf48
Binary files /dev/null and b/docs/source/images/au_map.png differ
diff --git a/docs/source/images/f_00.png b/docs/source/images/f_00.png
new file mode 100644
index 0000000000000000000000000000000000000000..833dde09baa0ea5d077bb0738455d08b353f82eb
Binary files /dev/null and b/docs/source/images/f_00.png differ
diff --git a/docs/source/images/f_01.png b/docs/source/images/f_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..3389363f2dba32e600c468cc944589fd22ac026c
Binary files /dev/null and b/docs/source/images/f_01.png differ
diff --git a/docs/source/images/f_02.png b/docs/source/images/f_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..bfc23d984fdc59c2da1b241067030fa40e899182
Binary files /dev/null and b/docs/source/images/f_02.png differ
diff --git a/docs/source/images/f_03.png b/docs/source/images/f_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..bdfa5d2969750945a91cef49d01e1eb79439147b
Binary files /dev/null and b/docs/source/images/f_03.png differ
diff --git a/docs/source/images/f_04.png b/docs/source/images/f_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..5bc47381caed3f0d8e8484eddd9aab711c572db6
Binary files /dev/null and b/docs/source/images/f_04.png differ
diff --git a/docs/source/images/f_05.png b/docs/source/images/f_05.png
new file mode 100644
index 0000000000000000000000000000000000000000..f5e01d3834cf7fb9137b09e37174ee0b7e61bbc7
Binary files /dev/null and b/docs/source/images/f_05.png differ
diff --git a/docs/source/images/hp_d.png b/docs/source/images/hp_d.png
new file mode 100644
index 0000000000000000000000000000000000000000..6a917412ec6ffee3fc0fd68a4c13edbd270057a3
Binary files /dev/null and b/docs/source/images/hp_d.png differ
diff --git a/docs/source/images/hp_m.png b/docs/source/images/hp_m.png
new file mode 100644
index 0000000000000000000000000000000000000000..0ef137f848afe2164aacae06e917f811a645622b
Binary files /dev/null and b/docs/source/images/hp_m.png differ
diff --git a/docs/source/images/hp_s.png b/docs/source/images/hp_s.png
new file mode 100644
index 0000000000000000000000000000000000000000..cd1ff14213e4dba261199907305f975778c1225b
Binary files /dev/null and b/docs/source/images/hp_s.png differ
diff --git a/docs/source/images/lens_d0.png b/docs/source/images/lens_d0.png
new file mode 100644
index 0000000000000000000000000000000000000000..533ba0b22589cc402d10e578a5a94b49d16d49e3
Binary files /dev/null and b/docs/source/images/lens_d0.png differ
diff --git a/docs/source/images/lens_m1.png b/docs/source/images/lens_m1.png
new file mode 100644
index 0000000000000000000000000000000000000000..2b0f579c900eb5c8167ab8b480940eca4606e2a0
Binary files /dev/null and b/docs/source/images/lens_m1.png differ
diff --git a/docs/source/images/lens_m2.png b/docs/source/images/lens_m2.png
new file mode 100644
index 0000000000000000000000000000000000000000..9f6e4f8859e15c8f5b27213323561620a77f4e0a
Binary files /dev/null and b/docs/source/images/lens_m2.png differ
diff --git a/docs/source/images/lens_s0.png b/docs/source/images/lens_s0.png
new file mode 100644
index 0000000000000000000000000000000000000000..a2281277d41489e38bbbd4e58574a7afd7ec2391
Binary files /dev/null and b/docs/source/images/lens_s0.png differ
diff --git a/docs/source/images/moon_d.png b/docs/source/images/moon_d.png
new file mode 100644
index 0000000000000000000000000000000000000000..621b6a8e99c3ac83ca4854e0851ca58ece65bcb2
Binary files /dev/null and b/docs/source/images/moon_d.png differ
diff --git a/docs/source/images/moon_kernel.png b/docs/source/images/moon_kernel.png
new file mode 100644
index 0000000000000000000000000000000000000000..b2104608c442fcbd8a8e15c361fd1ff6608dd787
Binary files /dev/null and b/docs/source/images/moon_kernel.png differ
diff --git a/docs/source/images/moon_m.png b/docs/source/images/moon_m.png
new file mode 100644
index 0000000000000000000000000000000000000000..89654fe47f9018131845b560398f6629ab55ccdd
Binary files /dev/null and b/docs/source/images/moon_m.png differ
diff --git a/docs/source/images/moon_mask.png b/docs/source/images/moon_mask.png
new file mode 100644
index 0000000000000000000000000000000000000000..1afbf06f7aa4b68ddd7d716130a7bb231d55e48d
Binary files /dev/null and b/docs/source/images/moon_mask.png differ
diff --git a/docs/source/images/moon_s.png b/docs/source/images/moon_s.png
new file mode 100644
index 0000000000000000000000000000000000000000..b36c1ff934d5898a9833de3d24f38622964ec52b
Binary files /dev/null and b/docs/source/images/moon_s.png differ
diff --git a/docs/source/images/moon_sigma.png b/docs/source/images/moon_sigma.png
new file mode 100644
index 0000000000000000000000000000000000000000..824c8157ba3d52cc09bb34dc672b42207c5abdd1
Binary files /dev/null and b/docs/source/images/moon_sigma.png differ
diff --git a/docs/source/images/radio_CLEAN.png b/docs/source/images/radio_CLEAN.png
new file mode 100644
index 0000000000000000000000000000000000000000..dc05d908caec0e43c3a305d3b9c06da16a81f9bb
Binary files /dev/null and b/docs/source/images/radio_CLEAN.png differ
diff --git a/docs/source/images/radio_RESOLVE.png b/docs/source/images/radio_RESOLVE.png
new file mode 100644
index 0000000000000000000000000000000000000000..b0a84022e74b42e3d5838687c66f41b53e904db1
Binary files /dev/null and b/docs/source/images/radio_RESOLVE.png differ
diff --git a/docs/source/images/radio_signal.png b/docs/source/images/radio_signal.png
new file mode 100644
index 0000000000000000000000000000000000000000..f51176c5aeb9f1ecd2933cd775ec6a78485721b1
Binary files /dev/null and b/docs/source/images/radio_signal.png differ
diff --git a/docs/source/images/rg1_d.png b/docs/source/images/rg1_d.png
new file mode 100644
index 0000000000000000000000000000000000000000..ab531ac499c6c0d904e0197fb79e98397a03b3b2
Binary files /dev/null and b/docs/source/images/rg1_d.png differ
diff --git a/docs/source/images/rg1_d_gap.png b/docs/source/images/rg1_d_gap.png
new file mode 100644
index 0000000000000000000000000000000000000000..92467cae0c450444c46ef9d2b7b7841934d6da71
Binary files /dev/null and b/docs/source/images/rg1_d_gap.png differ
diff --git a/docs/source/images/rg1_m.png b/docs/source/images/rg1_m.png
new file mode 100644
index 0000000000000000000000000000000000000000..daed9a990a2799f73a1dc19e8d77120cf15f1f30
Binary files /dev/null and b/docs/source/images/rg1_m.png differ
diff --git a/docs/source/images/rg1_m_err_.png b/docs/source/images/rg1_m_err_.png
new file mode 100644
index 0000000000000000000000000000000000000000..1f64cf31650d04cb07e5641d684663eef8fbc245
Binary files /dev/null and b/docs/source/images/rg1_m_err_.png differ
diff --git a/docs/source/images/rg1_m_gap_err_.png b/docs/source/images/rg1_m_gap_err_.png
new file mode 100644
index 0000000000000000000000000000000000000000..03afdd42f51c483688bd2df9c67e01a6213584c0
Binary files /dev/null and b/docs/source/images/rg1_m_gap_err_.png differ
diff --git a/docs/source/images/rg1_s.png b/docs/source/images/rg1_s.png
new file mode 100644
index 0000000000000000000000000000000000000000..a7d03d1a028ef0177994ae14d21c68e7befde563
Binary files /dev/null and b/docs/source/images/rg1_s.png differ
diff --git a/docs/source/images/rg2_d_pm.png b/docs/source/images/rg2_d_pm.png
new file mode 100644
index 0000000000000000000000000000000000000000..8a70dacd9a6374ceecb4396a6d283f195d66d309
Binary files /dev/null and b/docs/source/images/rg2_d_pm.png differ
diff --git a/docs/source/images/rg2_m_pm.png b/docs/source/images/rg2_m_pm.png
new file mode 100644
index 0000000000000000000000000000000000000000..ce7c9210baa16d674062dd44ae2f955c9b7e0d99
Binary files /dev/null and b/docs/source/images/rg2_m_pm.png differ
diff --git a/docs/source/images/rg2_s_pm.png b/docs/source/images/rg2_s_pm.png
new file mode 100644
index 0000000000000000000000000000000000000000..2dcae29c842f1a0c64948633d1096aa7a4856867
Binary files /dev/null and b/docs/source/images/rg2_s_pm.png differ
diff --git a/docs/source/images/t_00.png b/docs/source/images/t_00.png
new file mode 100644
index 0000000000000000000000000000000000000000..07880a2717c20ebe5140a4f0e015e1433833db6e
Binary files /dev/null and b/docs/source/images/t_00.png differ
diff --git a/docs/source/images/t_01.png b/docs/source/images/t_01.png
new file mode 100644
index 0000000000000000000000000000000000000000..83d4d7a24931badd1b7f480e437e9c21098d19cf
Binary files /dev/null and b/docs/source/images/t_01.png differ
diff --git a/docs/source/images/t_02.png b/docs/source/images/t_02.png
new file mode 100644
index 0000000000000000000000000000000000000000..9e6d288009a50d540bdd12d01e8563d45a1e32e9
Binary files /dev/null and b/docs/source/images/t_02.png differ
diff --git a/docs/source/images/t_03.png b/docs/source/images/t_03.png
new file mode 100644
index 0000000000000000000000000000000000000000..b72fcb6328f779b0f7ce8d4902b7f3111bfcc104
Binary files /dev/null and b/docs/source/images/t_03.png differ
diff --git a/docs/source/images/t_04.png b/docs/source/images/t_04.png
new file mode 100644
index 0000000000000000000000000000000000000000..f263b990c0202e357e4aedcf0e974e59a16ab821
Binary files /dev/null and b/docs/source/images/t_04.png differ
diff --git a/docs/source/images/t_05.png b/docs/source/images/t_05.png
new file mode 100644
index 0000000000000000000000000000000000000000..f60390555b938aa7ab4ea61c26532eeb4fcf7103
Binary files /dev/null and b/docs/source/images/t_05.png differ
diff --git a/docs/source/images/t_06.png b/docs/source/images/t_06.png
new file mode 100644
index 0000000000000000000000000000000000000000..5e28968932981368359a502a92a490a8135ee3fc
Binary files /dev/null and b/docs/source/images/t_06.png differ
diff --git a/docs/source/index.rst b/docs/source/index.rst
index bd9d79b8e3f97ff19b8c4a2da4555c1b258cffaa..f65d7a9e109810300c858e3bb9c2a8e8355bfeaf 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -23,9 +23,12 @@ Contents
 .. toctree::
    :maxdepth: 1
    
-   spaces
-   operator
-   minimization
+   gallery
+   ift
+   start
+   spaces/spaces
+   operators/operator
+   minimizer/minimization
    field
 
 
@@ -33,6 +36,6 @@ Indices and tables
 ..................
 
 * :ref:`genindex`
-* :ref:`modindex`
+* :any:`Module Index <mod/modules>`
 * :ref:`search`
 
diff --git a/docs/source/invertible_operator_mixin.rst b/docs/source/invertible_operator_mixin.rst
deleted file mode 100644
index 40d7392bd30f5d395a323d1bb8a131af384dcc73..0000000000000000000000000000000000000000
--- a/docs/source/invertible_operator_mixin.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``InvertibleOperatorMixin`` class -- ...
-............................................
-
-.. autoclass:: InvertibleOperatorMixin
-    :show-inheritance:
-    :members:
diff --git a/docs/source/conjugate_gradient.rst b/docs/source/minimizer/conjugate_gradient.rst
similarity index 100%
rename from docs/source/conjugate_gradient.rst
rename to docs/source/minimizer/conjugate_gradient.rst
diff --git a/docs/source/descent_minimizer.rst b/docs/source/minimizer/descent_minimizer.rst
similarity index 100%
rename from docs/source/descent_minimizer.rst
rename to docs/source/minimizer/descent_minimizer.rst
diff --git a/docs/source/minimization.rst b/docs/source/minimizer/minimization.rst
similarity index 100%
rename from docs/source/minimization.rst
rename to docs/source/minimizer/minimization.rst
diff --git a/docs/source/composed_operator.rst b/docs/source/operators/composed_operator.rst
similarity index 64%
rename from docs/source/composed_operator.rst
rename to docs/source/operators/composed_operator.rst
index db744c386f862aec6c02f9b9c5c0b12a50aac3c4..fea911d71c686009cb70a51146e17f367741f473 100644
--- a/docs/source/composed_operator.rst
+++ b/docs/source/operators/composed_operator.rst
@@ -1,7 +1,7 @@
 .. currentmodule:: nifty
 
-The ``ComposedOperator`` class -- A possibility to combine several operators
-............................................................................
+The ``ComposedOperator`` class -- NIFTY class for composed operators
+....................................................................
 
 .. autoclass:: ComposedOperator
     :show-inheritance:
diff --git a/docs/source/operators/diagonal_operator.rst b/docs/source/operators/diagonal_operator.rst
new file mode 100644
index 0000000000000000000000000000000000000000..f06eaf777350fa07526f1c813385ea2e27c43179
--- /dev/null
+++ b/docs/source/operators/diagonal_operator.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``DiagonalOperator`` class -- NIFTY class for diagonal operators
+....................................................................
+
+.. autoclass:: DiagonalOperator
+    :show-inheritance:
+    :members:
diff --git a/docs/source/endomorphic_operator.rst b/docs/source/operators/endomorphic_operator.rst
similarity index 80%
rename from docs/source/endomorphic_operator.rst
rename to docs/source/operators/endomorphic_operator.rst
index eb03cb726fbfb60006e2cd23c341934c549c6aee..facb826d54c6019ccbda28df5c0261ea4e898854 100644
--- a/docs/source/endomorphic_operator.rst
+++ b/docs/source/operators/endomorphic_operator.rst
@@ -3,6 +3,6 @@
 The ``EndomorphicOperator`` class -- ...
 ........................................
 
-.. autoclass:: EnomorphicOperator
+.. autoclass:: EndomorphicOperator
     :show-inheritance:
     :members:
diff --git a/docs/source/operators/fft_operator.rst b/docs/source/operators/fft_operator.rst
new file mode 100644
index 0000000000000000000000000000000000000000..92ac3b30f5a6b66e818addc5f281223e090daea2
--- /dev/null
+++ b/docs/source/operators/fft_operator.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``FFTOperator`` class -- Transforms between a pair of position and harmonic domains
+.......................................................................................
+
+.. autoclass:: FFTOperator
+    :show-inheritance:
+    :members:
diff --git a/docs/source/operators/invertible_operator_mixin.rst b/docs/source/operators/invertible_operator_mixin.rst
new file mode 100644
index 0000000000000000000000000000000000000000..3823b68ce00798306e4aa5dacb71200c4e009521
--- /dev/null
+++ b/docs/source/operators/invertible_operator_mixin.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``InvertibleOperatorMixin`` class -- Mixin class to invert implicit defined operators
+.........................................................................................
+
+.. autoclass:: InvertibleOperatorMixin
+    :show-inheritance:
+    :members:
diff --git a/docs/source/operator.rst b/docs/source/operators/operator.rst
similarity index 100%
rename from docs/source/operator.rst
rename to docs/source/operators/operator.rst
diff --git a/docs/source/operators/projection_operator.rst b/docs/source/operators/projection_operator.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d6985940120b14f74ea23310fa1b6fb44387fd27
--- /dev/null
+++ b/docs/source/operators/projection_operator.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``ProjectionOperator`` class -- NIFTY class for projection operators
+........................................................................
+
+.. autoclass:: ProjectionOperator
+    :show-inheritance:
+    :members:
diff --git a/docs/source/operators/propagator_operator.rst b/docs/source/operators/propagator_operator.rst
new file mode 100644
index 0000000000000000000000000000000000000000..ca468836efa43cd2628e4136f2f7dc2318bda7f8
--- /dev/null
+++ b/docs/source/operators/propagator_operator.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``PropagatorOperator`` class -- NIFTY Propagator Operator D
+...............................................................
+
+.. autoclass:: PropagatorOperator
+    :show-inheritance:
+    :members:
diff --git a/docs/source/response_operator.rst b/docs/source/operators/response_operator.rst
similarity index 68%
rename from docs/source/response_operator.rst
rename to docs/source/operators/response_operator.rst
index 4cfc9fed1fa7b80ccef2befdcc8bec316cdabd50..d55c178f3b93ea1a16b00af43dab4091fc4ffacc 100644
--- a/docs/source/response_operator.rst
+++ b/docs/source/operators/response_operator.rst
@@ -1,7 +1,7 @@
 .. currentmodule:: nifty
 
-The ``ResponseOperator`` class -- A possible response implementation
-....................................................................
+The ``ResponseOperator`` class -- NIFTy ResponseOperator (example)
+..................................................................
 
 .. autoclass:: ResponseOperator
     :show-inheritance:
diff --git a/docs/source/operators/smoothing_operator.rst b/docs/source/operators/smoothing_operator.rst
new file mode 100644
index 0000000000000000000000000000000000000000..68057813da81f8209b37dbc85404ea379a650ae8
--- /dev/null
+++ b/docs/source/operators/smoothing_operator.rst
@@ -0,0 +1,8 @@
+.. currentmodule:: nifty
+
+The ``SmoothingOperator`` class -- NIFTY class for smoothing operators
+......................................................................
+
+.. autoclass:: SmoothingOperator
+    :show-inheritance:
+    :members:
diff --git a/docs/source/projection_operator.rst b/docs/source/projection_operator.rst
deleted file mode 100644
index 761f4373e601bb6f8b1982bba86aeb7a95c3d332..0000000000000000000000000000000000000000
--- a/docs/source/projection_operator.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``ProjectionOperator`` class -- ...
-.......................................
-
-.. autoclass:: ProjectionOperator
-    :show-inheritance:
-    :members:
diff --git a/docs/source/propagator_operator.rst b/docs/source/propagator_operator.rst
deleted file mode 100644
index 28fcbcf9a20cfcb933ea366a2dbd5e3125e1fd18..0000000000000000000000000000000000000000
--- a/docs/source/propagator_operator.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``PropagatorOperator`` class -- ...
-.......................................
-
-.. autoclass:: PropagatorOperator
-    :show-inheritance:
-    :members:
diff --git a/docs/source/rg_transforms.rst b/docs/source/rg_transforms.rst
deleted file mode 100644
index 54b3a0ce535cbc21658f2b37549d18ce7c84caf9..0000000000000000000000000000000000000000
--- a/docs/source/rg_transforms.rst
+++ /dev/null
@@ -1,9 +0,0 @@
-.. currentmodule:: nifty
-
-The ``Transform`` class -- A transformation routine
-...................................................
-
-
-.. autoclass:: Transform
-    :show-inheritance:
-    :members:
diff --git a/docs/source/smoothing_operator.rst b/docs/source/smoothing_operator.rst
deleted file mode 100644
index 39f8a403038b0042301c41f16ff2d434b400eb47..0000000000000000000000000000000000000000
--- a/docs/source/smoothing_operator.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-.. currentmodule:: nifty
-
-The ``SmoothingOperator`` class -- Smoothing fields
-...................................................
-
-.. autoclass:: SmoothingOperator
-    :show-inheritance:
-    :members:
diff --git a/docs/source/gl_space.rst b/docs/source/spaces/gl_space.rst
similarity index 100%
rename from docs/source/gl_space.rst
rename to docs/source/spaces/gl_space.rst
diff --git a/docs/source/gllm_transformation.rst b/docs/source/spaces/gllm_transformation.rst
similarity index 100%
rename from docs/source/gllm_transformation.rst
rename to docs/source/spaces/gllm_transformation.rst
diff --git a/docs/source/hp_space.rst b/docs/source/spaces/hp_space.rst
similarity index 100%
rename from docs/source/hp_space.rst
rename to docs/source/spaces/hp_space.rst
diff --git a/docs/source/hplm_transformation.rst b/docs/source/spaces/hplm_transformation.rst
similarity index 100%
rename from docs/source/hplm_transformation.rst
rename to docs/source/spaces/hplm_transformation.rst
diff --git a/docs/source/lm_space.rst b/docs/source/spaces/lm_space.rst
similarity index 100%
rename from docs/source/lm_space.rst
rename to docs/source/spaces/lm_space.rst
diff --git a/docs/source/lmgl_transformation.rst b/docs/source/spaces/lmgl_transformation.rst
similarity index 100%
rename from docs/source/lmgl_transformation.rst
rename to docs/source/spaces/lmgl_transformation.rst
diff --git a/docs/source/lmhp_transformation.rst b/docs/source/spaces/lmhp_transformation.rst
similarity index 100%
rename from docs/source/lmhp_transformation.rst
rename to docs/source/spaces/lmhp_transformation.rst
diff --git a/docs/source/power_space.rst b/docs/source/spaces/power_space.rst
similarity index 65%
rename from docs/source/power_space.rst
rename to docs/source/spaces/power_space.rst
index 1ed47f4a9d07b14b50378c909579e503753c0802..589cd011748e36858d2b9006e34dc3d553ebd7b2 100644
--- a/docs/source/power_space.rst
+++ b/docs/source/spaces/power_space.rst
@@ -1,7 +1,7 @@
 .. currentmodule:: nifty
 
-The ``PowerSpace`` class --  The natural space underlying power-spectra
-.......................................................................
+The ``PowerSpace`` class --  NIFTY class for spaces of power spectra
+....................................................................
 
 .. autoclass:: PowerSpace
     :show-inheritance:
diff --git a/docs/source/rg_space.rst b/docs/source/spaces/rg_space.rst
similarity index 100%
rename from docs/source/rg_space.rst
rename to docs/source/spaces/rg_space.rst
diff --git a/docs/source/rgrg_transformation.rst b/docs/source/spaces/rgrg_transformation.rst
similarity index 100%
rename from docs/source/rgrg_transformation.rst
rename to docs/source/spaces/rgrg_transformation.rst
diff --git a/docs/source/spaces.rst b/docs/source/spaces/spaces.rst
similarity index 98%
rename from docs/source/spaces.rst
rename to docs/source/spaces/spaces.rst
index 43c0c61bdc3929a478269204b40bb89a7bc11ad7..c918b316c012f56616af1815ca7108e2a7967961 100644
--- a/docs/source/spaces.rst
+++ b/docs/source/spaces/spaces.rst
@@ -19,6 +19,7 @@ Next to the generic :py:class:`Space` class, NIFTY has implementations of five s
     gl_space
     lm_space
     power_space
+    transformations
 
 .. currentmodule:: nifty
 
diff --git a/docs/source/transformations.rst b/docs/source/spaces/transformations.rst
similarity index 75%
rename from docs/source/transformations.rst
rename to docs/source/spaces/transformations.rst
index 153a7689d448fd061d8a4b20e7ea4c465007fd5a..47450a6ac64e291f6c4d8b343df1fb69f0628aec 100644
--- a/docs/source/transformations.rst
+++ b/docs/source/spaces/transformations.rst
@@ -1,11 +1,10 @@
 Transformations
 ---------------
-NIFTY provides transformations
+NIFTY helper classes for transformations betwenn NIFTY spaces.
 
 .. toctree:: 
     :maxdepth: 1
     
-    rg_transforms
     rgrg_transformation
     gllm_transformation
     hplm_transformation
diff --git a/docs/source/start.rst b/docs/source/start.rst
new file mode 100644
index 0000000000000000000000000000000000000000..100947055000e43511ac79aaafa91a1ea3f0ef04
--- /dev/null
+++ b/docs/source/start.rst
@@ -0,0 +1,784 @@
+.. _informal_label:
+
+First steps -- An informal introduction
+=======================================
+
+New to Python?
+--------------
+
+If you are new to Python, it is recommended that you familiarize yourself with Python before proceeding with the NIFTY tutorial. Start with `an informal introduction to Python <http://docs.python.org/2/tutorial/introduction.html>`_.
+
+NIFTY Tutorial
+--------------
+
+.. currentmodule:: nifty
+
+.. automodule:: nifty
+
+NIFTY enables the programming of grid and resolution independent algorithms. In particular for signal inference algorithms, where a continuous signal field is to be recovered, this freedom is desirable. This is achieved with an object-oriented infrastructure that comprises, among others, abstract classes for :ref:`spaces <spaces>`, :ref:`fields <fields>`, and :ref:`operators <operators>`. All those are covered in this tutorial.
+
+You should be able to import NIFTY like this after a successfull `installation <install.html>`_.
+
+>>> from nifty import *
+
+.. note:: During the import of NIFTY, other (common) modules get imported as well: ``import os``, ``import numpy as np``, and ``import pylab as pl``.
+
+For a version statement and further settings you can access the global ``about`` variable.
+
+>>> about
+nifty version 0.4.0
+
+.. _spaces:
+
+Spaces
+......
+
+The very foundation of NIFTY's framework are the spaces, of which there are:
+
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+| Space subclass           | Corresponding grid                                                | Conjugate space class                        |
++==========================+===================================================================+==============================================+
+| :py:class:`Space`        | the base space class                                              | (none)                                       |
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+| :py:class:`RGSpace`      | *n*-dimensional regular Euclidean grid over :math:`\mathcal{T}^n` | :py:class:`RGSpace`                          |
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+| :py:class:`LMSpace`      | spherical harmonics                                               | :py:class:`GLSpace` or :py:class:`HPSpace`   |
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+| :py:class:`GLSpace`      | Gauss-Legendre grid on the :math:`\mathcal{S}^2` sphere           | :py:class:`LMSpace`                          |
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+| :py:class:`HPSpace`      | HEALPix grid on the :math:`\mathcal{S}^2` sphere                  | :py:class:`LMSpace`                          |
++--------------------------+-------------------------------------------------------------------+----------------------------------------------+
+
+
+A space instance can be created by providing all required parameters. The printout of a space informs you about all the parameters set.
+
+
+>>> space = RGSpace(10, naxes=2) # two-dimensional regular 10 x 10 grid of physical size 1 in position space
+>>> print(space)
+nifty.RGSpace instance
+- num        = [10, 10]
+- naxes      = 2
+- hermitian  = True
+- purelyreal = True
+- zerocenter = [True, True]
+- dist       = [0.1, 0.1]
+- fourier    = False
+
+>>> space = HPSpace(32) # HEALPix grid with nside = 32
+>>> print(space)
+nifty.HPSpace instance
+- nside = 32
+
+.. note:: All instances of space classes can be compared by ``==`` or ``!=`` (or ``<>``); instances of the same class allow additionally for comparison by ``<``, ``<=``, ``>=``, or ``>``.
+
+The dimensionality of a space and its number of degrees of freedom are fundamental characteristics of a space, and can be accessed by :py:meth:`space.dim` and :py:meth:`space.dof`. Fields, which are discussed later, are defined on spaces and the array in which the field values are stored will be of the "split" dimensionality of this space.
+
+
+>>> space = RGSpace(10, naxes=2)
+>>> space.dim(split=False) # 100 grid points
+100
+>>> space.dim(split=True) # structured 10 x 10
+array([10, 10])
+
+>>> space = HPSpace(32)
+>>> space.dim() # == 12 * nside ** 2
+12288
+
+Nested spaces
+,,,,,,,,,,,,,
+
+The :py:class:`nested_space` class is designed to handle arbitrary product spaces. This is useful for outer products of fields and for grouping fields.
+
+>>> space = nested_space([point_space(8), rg_space(10, naxes=2), hp_space(32)])
+>>> print(space)
+nifty.nested_space instance
+- nest = [<nifty.point_space>, <nifty.rg_space>, <nifty.hp_space>]
+>>> space.dim(split=True)
+array([    8,    10,    10, 12288])
+
+Conjugate spaces
+,,,,,,,,,,,,,,,,
+
+Spatial symmetries of a system can be exploited by corresponding coordinate transformations. Often, transformations from one basis to its harmonic counterpart can greatly reduce the computational complexity of an algorithm. The "harmonic" (or "conjugate") basis is defined by the eigenbasis of the Laplace operator. This conjugation of bases is implemented in NIFTY by distinguishing conjugate space classes.
+
+Since the basis transformation from one basis to itself is trivial, the evaluation of :py:meth:`space.check_codomain` to itself always returns ``True``.
+
+
+A valid conjugate space is obtained by :py:meth:`space.get_codomain`. For example, the conjugate basis of a flat position space is the Fourier basis.
+
+>>> x_space = RGSpace(10, dist=0.2) # position space housing purely real field values
+>>> print(x_space)
+nifty.RGSpace instance
+- num        = [10]
+- naxes      = 1
+- hermitian  = True
+- purelyreal = True
+- zerocenter = [True]
+- dist       = [0.2]
+- fourier    = False
+>>> k_space = x_space.get_codomain() # Fourier space housing complex field values
+>>> print(k_space)
+nifty.rg_space instance
+- num        = [10]
+- naxes      = 1
+- hermitian  = True
+- purelyreal = False
+- zerocenter = [True]
+- dist       = [0.5]
+- fourier    = True
+>>> x_space.check_codomain(k_space)
+True
+>>> x_space == k_space.get_codomain()
+True
+>>> x_space.dof() == k_space.dof() == 10 # due to hermitian symmetry
+True
+
+For spherical position spaces a spherical harmonics transformation yields the conjugate basis.
+
+>>> k_space = lm_space(95) # spherical harmonics with lmax = 95
+>>> print(k_space)
+nifty.lm_space instance
+- lmax     = 95
+- mmax     = 95
+- datatype = numpy.complex128
+>>> print(k_space.get_codomain(coname="gl")) # corresponding Gauss-Legendre grid
+nifty.gl_space instance
+- nlat     = 96
+- nlon     = 191
+- datatype = numpy.float64
+>>> print(k_space.get_codomain(coname="hp")) # corresponding HEALPix grid
+nifty.hp_space instance
+- nside = 32
+
+The physical volume of a conceptually continuous space is returned by :py:meth:`space.get_meta_volume` (per pixel or in total). This is not useful for discrete spaces, such as the :py:class:`point_space` or the :py:class:`lm_space`.
+
+>>> space = point_space(8)
+>>> space.discrete
+True
+
+>>> x_space = rg_space(10, dist=0.2)
+>>> x_space.get_meta_volume(total=True) # == 10*0.2
+2.0
+>>> k_space = x_space.get_codomain()
+>>> k_space.dist() # == 1/(10*0.2)
+array([ 0.5])
+>>> k_space.get_meta_volume(total=True) # == 10*0.5
+5.0
+
+>>> space = hp_space(32)
+>>> space.get_meta_volume(total=True) # == 4*pi
+12.566370614359172
+
+.. _power_index:
+
+`(Read more about spaces.) <space.html>`_
+
+Power indexing **(advanced)**
+,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+
+The basis transformations mapping positions to wave vectors, :math:`\vec{x} = (x_1,x_2,\dots) \mapsto \vec{k} = (k_1,k_2,\dots)`, contain several computationally interesting quantities that are summarized in the following. Those quantities, called "power indices", can be accessed in different ways: (a) by the :py:attr:`space.power_indices` dictionary after the :py:meth:`space.set_power_indices` method has been evaluated once, (b) automatically by forwarding the respective space instance, or (c) by handling the output of the :py:meth:`space.get_power_indices` method; the latter is not recommended. (Only spaces representing the conjugate basis of a position space can evaluate those instance methods without raising an exception.)
+
+For each grid point position :math:`\vec{k}` one can compute its abolute :math:`k = \big| \vec{k} \big|`. All the appearing absolute values (*irreducible spectral indices*) and their frequency of appearence (*number of degrees of freedom*) are represented as ``kindex`` and ``rho`` for each *irreducible spectral band*. In turn, one can assign each band an integer indicating the position in the ``kindex`` array and thus fill an array of the full size of the space that is denoted ``pindex``. (Indexing with the *unindexing* list ``pundex`` undoes the indexing with the *indexing* array ``pindex``.)
+
+How does this look in practice? Below, a regular :math:`4 \times 4` grid is initialized with unit distance. The ``pindex`` shows :math:`4 \times 4` array with ``0`` at its "center" :math:`(k=0)`. The four nearest neighbors :math:`(k=1)` are right next to it and four second next neighbors :math:`(k=\sqrt{2})` touch the edges of the "center". Accordingly, the ``kindex`` reads :math:`(0,1,\sqrt{2},\dots)` and the ``rho`` :math:`(1,4,4,\dots)`.
+
+>>> k_space = rg_space([4, 4], dist=1, fourier=True) # Fourier space
+>>> k_space.set_power_indices() # power indices are computed
+>>> k_space.power_indices["pindex"]
+array([[5, 4, 3, 4],
+       [4, 2, 1, 2],
+       [3, 1, 0, 1],
+       [4, 2, 1, 2]])
+>>> k_space.power_indices["kindex"]
+array([ 0.        ,  1.        ,  1.41421356,  2.        ,  2.23606798,  2.82842712])
+>>> k_space.power_indices["rho"]
+array([1, 4, 4, 2, 4, 1])
+
+.. note:: The :py:meth:`rg_space.set_power_indices` method allows for a (re)binning of the spectral indices which is useful for high dimensional settings.
+
+Since the :py:class:`lm_space` is discrete, the angular quantum number :math:`\ell` (instead of some :math:`\big| \vec{\ell} \big|`) labels the indices of which there are :math:`2\ell+1`.
+
+>>> k_space = lm_space(3) # spherical harmonics with lmax = 3
+>>> k_space.set_power_indices() # power indices are computed
+>>> k_space.power_indices["pindex"]
+array([0, 1, 2, 3, 1, 2, 3, 2, 3, 3])
+>>> k_space.power_indices["kindex"]
+array([0, 1, 2, 3])
+>>> k_space.power_indices["rho"]
+array([1, 3, 5, 7])
+
+One advantage of basis conjugation is that the covariance of a :ref:`Gaussian random field <Garfields>` that is statistically homogeneous in position space becomes diagonal in the harmonic basis. This diagonal is then called power spectrum. Therefore, operations that involve power spectra or spectral bands in general require "power indexing".
+
+.. _fields:
+
+Fields
+......
+
+The actual algorithmic work is done on the level of fields. In order to initialize a field instance, a ``domain``, the space in which the field is defined, needs to be stated. The field values ``val`` are all set to zero, if not given, and the ``target``, the default conjugate space, is obtained by default using ``domain.get_codomain()``. The printout of a field informs you about all those three properties of a field.
+
+>>> f = field(point_space(8)) # field of 8 points all carrying 0.0
+>>> print(f)
+nifty.field instance
+- domain      = <nifty.point_space>
+- val         = [...]
+  - min.,max. = [0.0, 0.0]
+  - med.,mean = [0.0, 0.0]
+- target      = <nifty.point_space>
+
+Any field can be initialized with random numbers by specifying the `random` keyword and further constraints.
+
+>>> f = field(rg_space([4, 4]), random="gau") # field of 4 x 4 grid points filled with Gaussian random numbers
+>>> f.val
+array([[-0.48543864, -1.56595855,  0.38683479,  1.08021123],
+       [ 0.26197572, -0.2245229 ,  1.61690895, -1.28385473],
+       [ 3.04056877, -0.48699409,  0.44017361,  0.03244873],
+       [-0.88216182, -0.11849966,  0.51377166, -0.37855755]])
+>>> f.domain.fourier
+False
+>>> f.target.fourier
+True
+
+More informative than the printout of a field is of course a plot obtained by :py:meth:`field.plot`.
+
+>>> space = hp_space(32)
+>>> f = field(space, val=np.arange(space.dim()))
+>>> f.plot() # a pylab figure will open.
+
+.. image:: images/t_00.png
+    :width: 33 %
+
+Basic operations
+,,,,,,,,,,,,,,,,
+
+The field values are stored as a `numpy.ndarray <http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html>`_ and can be accessed by the ``val`` property or by indexing the field with ``[:]``.
+
+>>> f = field(point_space(8), val=3) # field of 8 points all carrying 3.0
+>>> f.val
+array([ 3.,  3.,  3.,  3.,  3.,  3.,  3.,  3.])
+>>> f[:] # standard numpy indexing
+array([ 3.,  3.,  3.,  3.,  3.,  3.,  3.,  3.])
+
+Standard operations (``+``, ``-``, ``*``, ``/``, ...) are available for fields of the same domain.
+
+>>> f = field(rg_space(4), val=3)
+>>> u = field(rg_space(4), val=[0, 1, 2, 3])
+>>> f + u
+<nifty.field>
+>>> (f ** u)[2] # == 3.0 ** 2.0
+9.0
+
+.. warning:: The Python interpreter reads from left to right. As a consequence, the compatibility to NumPy arrays introduces the following ambiguity, ``field + array == field`` but ``array + field == array``.
+
+Trigonometric, exponential, and logarithmic functions are also valid
+`field operations <field_operations.html>`_. They are applied for each
+field value individually.
+
+>>> f = field(hp_space(32), random="uni", vmin=-2, vmax=8)
+>>> exp(f).plot(vmin=exp(-2), vmax=exp(8)) # a pylab figure will open.
+
+.. image:: images/t_01.png
+    :width: 33 %
+
+Inner & outer products
+,,,,,,,,,,,,,,,,,,,,,,
+
+This section deals with the inner and outer products provided by :py:meth:`field.dot`, :py:meth:`field.pseudo_dot`, and :py:meth:`field.tensor_dot`. The crucial point is that those products involove position integrals over continuous position spaces that are approximated by weighted sums.
+
+The scalar product implemented in :py:meth:`field.dot` of two fields :math:`a` and :math:`b` describes the following approximation, where :math:`V_q` is the volume of the :math:`q`\th pixel (:math:`\Omega` some position space, :math:`x \in \Omega` a position, :math:`q \in \{1,\dots,Q\}` a discrete index; :math:`^\dagger` denotes transposition and complex conjugation, and :math:`^*` only the latter).
+
+.. math::
+
+    a^\dagger b \quad\equiv\quad \int_{\Omega} \mathrm{d}x \; a^*(x) \; b(x) \quad\approx\quad \sum_{q=1}^Q V_q \; a^*_q \; b_q
+
+In discrete spaces the scalar product is the common inner Euclidean product.
+
+>>> f = field(point_space(8), val=3) # here Vpix = 1
+>>> f.domain.discrete
+True
+>>> f.dot(f) # == 8 * (3.0 * 3.0)
+72.0
+
+In conceptually continuous spaces the scalar product respects the physical size of the space.
+
+>>> f = field(rg_space(4, dist=0.5), val=3) # here Vpix = 0.5
+>>> f.dot(f) # == 4 * 0.5 * (3.0 * 3.0)
+18.0
+
+>>> f = field(hp_space(32), val=0.25) # here Vpix = (4 * pi / Npix) and Npix = 12 * nside ** 2
+>>> f.dot(1) # == Npix * Vpix * (0.25 * 1) == pi
+3.1415926535896275
+>>> f.dot(1/f) # == Npix * Vpix * (0.25 * 1 / 0.25) == 4 * pi
+12.56637061435851
+
+.. note:: The relation :math:`V_\mathrm{pix} = V_\mathrm{total} / N_\mathrm{pix}` is not true for all spaces.
+
+The tensor product implemented in :py:meth:`field.tensor_dot` of two fields :math:`a` and :math:`b` yields a field :math:`c` that is defined on the product space of the domains of :math:`a` and :math:`b`. The tensor product is not symmetric.
+
+.. math::
+
+    a \otimes b = c \quad\equiv\quad c(x,y) = a(x) \cdot b(y)
+
+
+>>> a = field(point_space(8), val=np.arange(1, 9))
+>>> b = field(rg_space(4), val=[0, 2, 4, 8])
+>>> c = a.tensor_dot(b)
+>>> print(c)
+nifty.field instance
+- domain      = <nifty.nested_space>
+- val         = [...]
+  - min.,max. = [0.0, 64.0]
+  - med.,mean = [11.0, 15.75]
+- target      = <nifty.nested_space>
+>>> print(c.domain)
+nifty.nested_space instance
+- nest = [<nifty.point_space>, <nifty.rg_space>]
+>>> c.val # shape (8, 4)
+array([[  0.,   2.,   4.,   8.],
+       [  0.,   4.,   8.,  16.],
+       [  0.,   6.,  12.,  24.],
+       [  0.,   8.,  16.,  32.],
+       [  0.,  10.,  20.,  40.],
+       [  0.,  12.,  24.,  48.],
+       [  0.,  14.,  28.,  56.],
+       [  0.,  16.,  32.,  64.]])
+>>> b.tensor_dot(a)[:] # shape (4, 8)
+array([[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.],
+       [  2.,   4.,   6.,   8.,  10.,  12.,  14.,  16.],
+       [  4.,   8.,  12.,  16.,  20.,  24.,  28.,  32.],
+       [  8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.]])
+
+The pseudo scalar product implemented in :py:meth:`field.pseudo_dot` performs a scalar product in some dimensions and returns a scalar or a field, which might be explained best with an example following the last one.
+
+>>> a = field(point_space(8), val=np.arange(1, 9))
+>>> b = field(rg_space(4), val=[0, 2, 4, 8]) # Vpix = 0.25
+>>> c = a.tensor_dot(b)
+>>> b.dot(b) # == 0.25 * (0.0 ** 2 + 2.0 ** 2 + 4.0 ** 2 + 8.0 ** 2)
+21.0
+>>> c.pseudo_dot(b)
+<nifty.field>
+>>> c.pseudo_dot(b).val # == a * (b.dot(b))
+array([  21.,   42.,   63.,   84.,  105.,  126.,  147.,  168.])
+>>> b.dot(1) # == 0.25 * (0.0 * 1 + 2.0 * 1 + 4.0 * 1 + 8.0 * 1)
+3.5
+>>> c.pseudo_dot(1).val # == a * (b.dot(1))
+array([  3.5,   7. ,  10.5,  14. ,  17.5,  21. ,  24.5,  28. ])
+>>> c.dot(c) # == (a.dot(a)) * (b.dot(b))
+4284.0
+>>> c.dot(1) # == (a.dot(1)) * (b.dot(1))
+126.0
+
+Transformations
+,,,,,,,,,,,,,,,
+
+In order to take the full advantage of basis transformations, each field instance can evaluate its :py:meth:`field.transform` method that applies a transformation from the field's domain to its target. The field's transform is again a field and therefore a valid input for field methods. In favor of user-friendliness, transformations between harmonic bases may be applied automatically if required.
+
+>>> x_space = hp_space(32)
+>>> k_space = x_space.get_codomain()
+>>> f = field(x_space, val=1)
+>>> f.transform()
+<nifty.field>
+>>> f.domain == x_space
+True
+>>> f.target == k_space
+True
+>>> f.transform().domain == k_space
+True
+>>> f.transform().target == x_space
+True
+>>> f.dot(f)
+12.56637061435851
+>>> f.dot(f.transform())
+12.56637061435851
+>>> f.transform().dot(f)
+12.56637061435917
+>>> f.transform().dot(f.transform())
+12.56637061435917
+
+.. note:: Transformations between conjugate spaces are performed automatically if required.
+
+A non-standard transformation can be applied by explicitly stating a target that is a valid codomain (valid by means of `check_codomain` returns ``True``).
+
+>>> k_space = lm_space(383)
+>>> x1_space = k_space.get_codomain(coname="hp") # corresponding HEALPix grid
+>>> x2_space = k_space.get_codomain(coname="gl") # corresponding Gauss-Legendre grid
+>>> f = field(x1_space, val=np.load("/SOMEWHERE/nifty/demos/demo_faraday_map.npy")) # field values with some easily visible features
+>>> from nifty.nifty_cmaps import *
+>>> f.plot(vmin=-4, vmax=4, cmap=ncmap.fm()) # a pylab figure will open.
+>>> f.transform(target=k_space).transform(target=x2_space).plot(vmin=-4, vmax=4, cmap=ncmap.fm()) # another pylab figure will open.
+
++----------------------------+----------------------------+
+| .. image:: images/f_00.png | .. image:: images/f_01.png |
+|     :width:  66 %          |     :width:  66 %          |
++----------------------------+----------------------------+
+
+More methods
+,,,,,,,,,,,,
+
+The :py:class:`field` class offers further methods, which are listed in the following table, at your disposal. Some of those have been discussed already.
+
++------------------------------+-----------------------------------------------------------------------------------------------+
+| Method name                  | Description                                                                                   |
++==============================+===============================================================================================+
+| :py:meth:`Field.cast_domain` | Alters the field's domain without altering the field values or the codomain.                  |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.conjugate`   | Complex conjugates the field values.                                                          |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.dim`         | Returns the dimensionality of the field.                                                      |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.dot`         | Applies the scalar product between two fields, returns a scalar.                              |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.hat`         | Translates the field into a diagonal operator.                                                |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.inverse_hat` | Translates the inverted field into a diagonal operator.                                       |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.norm`        | Returns the :math:`L^2`-norm of the field.                                                    |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.plot`        | Draws a figure illustrating the field.                                                        |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.power`       | Computes the power spectrum of the field values.                                              |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.pseudo_dot`  | Applies a scalar product between two fields on a certain subspace of a product space, returns |
+|                              | a scalar or a field, depending on the subspace.                                               |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.set_target`  | Alters the field's codomain without altering the domain or the field values.                  |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.set_val`     | Alters the field values without altering the domain or codomain.                              |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.smooth`      | Smoothes the field values in position space by convolution with a Gaussian kernel.            |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.tensor_dot`  | Applies a tensor product between two fields, returns a field defined in the product space.    |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.transform`   | Applies a transformation from the field's domain to some codomain.                            |
++------------------------------+-----------------------------------------------------------------------------------------------+
+| :py:meth:`Field.weight`      | Multiplies the field with the grid's volume factors (to a given power).                       |
++------------------------------+-----------------------------------------------------------------------------------------------+
+
+`(Read more about fields) <field.html>`_
+
+.. _Garfields:
+
+Gaussian random fields **(advanced)**
+,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+
+Since NIFTY shall assist in the development of inference, Gaussian random fields are discussed in this section. Gaussian random fields are fields drawn from a `multivariate normal distribution <http://en.wikipedia.org/wiki/Multivariate_normal_distribution>`_ that is characterized by its mean and covariance. (Here, only zero-mean Gaussians are discussed.)
+
+The covariance is a symmetric and positive-definite operator of the squared dimension of the field. It describes the correlation between the field values: In its matrix representation, the covariance element :math:`(ij)` would describe the correlation between the field values :math:`i` and :math:`j`. The simplest covariance would be :math:`\delta_{ij}` producing field values that are uncorrelated, so called "white noise".
+
+>>> x_space = rg_space(100, naxes=2)
+>>> f = field(x_space, random="gau")
+>>> f.plot(vmin=-1, vmax=1, cmap=pl.cm.gray) # a pylab figure will open.
+
+.. image:: images/t_02.png
+    :width: 25 %
+
+If the covariance is only determined by a diagonal, this diagonal is referred to as *power spectrum*. The covariance, :math:`\delta_{ij}`, of the above field is denoted flat or "white" power spectrum. One way to introduce a correlation between the field values is to draw a field from a power spectrum that (a) bundles harmonic modes that have the same ``kindex``, :math:`k = \big| \vec{k} \big|`, and (b) has some non-flat structure. (Confer :ref:`power indexing <power_index>`.)
+
+The obtained field is statistically homogeneous and isotropic, that means the correlation between two field values at positions :math:`\vec{x}` and :math:`\vec{y}` depends only on their physical distance :math:`\big| \vec{x} - \vec{y} \big| \propto 1/k`. The power spectrum is returned by :py:meth:`field.power`, and can also be plotted by :py:meth:`field.plot`.
+
+>>> x_space = rg_space(100, naxes=2)
+>>> k_space = x_space.get_codomain()
+>>> spectrum_a = (lambda k: (1 / (k + 1)) ** 5)  # falling spectrum
+>>> spectrum_b = (lambda k: (3 / (k + 10)) ** 5) # peaking spectrum
+>>> a = field(x_space, target=k_space, random="syn", spec=spectrum_a)
+>>> b = field(x_space, target=k_space, random="syn", spec=spectrum_b)
+>>> a.plot()
+>>> a.plot(power=True, other=(spectrum_a,spectrum_b), mono=False, vmin=1E-5, vmax=1E-1)
+>>> b.plot()
+>>> b.plot(power=True, other=(spectrum_a,spectrum_b), mono=False, vmin=1E-5, vmax=1E-1)
+
++----------------------------+----------------------------+
+| .. image:: images/t_03.png | .. image:: images/t_04.png |
+|     :width:  50 %          |     :width:  70 %          |
++----------------------------+----------------------------+
+| .. image:: images/t_05.png | .. image:: images/t_06.png |
+|     :width:  50 %          |     :width:  70 %          |
++----------------------------+----------------------------+
+
+Note the difference between the given "theoretical" power spectrum (orange, blue) and the power spectrum from the field realization (green). A good way to get a feeling for Gaussian random fields is by playing around. Experiment with the spectra! (Adapt or omit the parameters ``vmin`` and ``vmax``!)
+
+.. note:: At this point you can grasp the idea of grid and resolution independent code. Try to change the position space; for example, set ``x_space = hp_space(32)``!
+
+.. _operators:
+
+Operators
+.........
+
+The missing ingredient for manufacturing algorithms is the :py:class:`operator` class. All operators are initialized with a ``domain`` (optionally also with a ``target``) and various flags. The ``domain`` specifies the space in which input fields are defined (and the ``target`` specifies the space for output fields). According to the requirements of the operators, other parameters, like the general ``para`` property, need to be stated.
+
+The most basic operator is returned by :py:func:`identity` and represents an identity in the specified domain. This (rather useless) operator shall serve as an educational example for getting to know the :py:class:`operator` class, or more precisely the :py:class:`diagonal_operator` class.
+
+>>> space = point_space(8)
+>>> f = field(space, val=3)
+>>> I = identity(space) # == identity operator
+>>> print(I)
+<nifty.diagonal_operator>
+>>> I.domain == I.target
+True
+>>> I(f) # == I.times(f) due to operator.__call__
+<nifty.field>
+>>> I.adjoint_inverse_times(f) # there are more operator.*times methods
+<nifty.field>
+>>> I(f).val
+array([ 3.,  3.,  3.,  3.,  3.,  3.,  3.,  3.])
+>>> I(5) # input may be raised to be a field
+<nifty.field>
+>>> I(5).val
+array([ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.])
+
+If the domain is a continuous space, the appropriate weighting is done automatically in order to correctly approximate position integrals. However, this trivial task implies several discretization subtleties. (It matters if one absorbs the weights into the "bare" matrix entries.)
+
+.. math::
+
+    \mathrm{id}[a(x)] \quad\equiv\quad \int_{\Omega} \mathrm{d}y \; \underbrace{\delta(x-y)}_\mathrm{continuous} \; a(y)
+    \quad\approx\quad \sum_{q=1}^Q V_q \; \underbrace{\phantom{P} \frac{\delta_{pq}}{V_q} \phantom{P}}_\mathrm{"bare"} \; a_q
+    \quad=\quad \sum_{q=1}^Q \; \underbrace{\phantom{P} \delta_{pq} \phantom{P}}_\mathrm{"non-bare"} \; a_q
+
+>>> space = rg_space(10, dist=0.2) # Vpix = 0.2
+>>> I = identity(space) # == identity operator
+>>> I.diag()
+array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])
+>>> I.diag(bare=True)
+array([ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.])
+>>> I.inverse_diag(bare=True)
+array([ 5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.,  5.])
+>>> I.hat() # diagonal to field
+<nifty.field>
+>>> I.hat().val
+array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.])
+
+The involvement of volume weights can cause some confusion concerning the interpretation of the corresponding matrix elements. As shown above, the discretization of the continuous identity operator, which equals a :math:`\delta`-distribution, yields a weighted Kronecker-Delta. Nevertheless, the "bare" entries are the ones to be consulted for interpretation. -- Say a field is drawn from a zero-mean Gaussian with a covariance that equals the identity. The intuitive assumption that the field values have a variance of :math:`1` is not true; the variance is given by :math:`1/V_\mathrm{pix}`. Convince yourself by drawing a field with ``I.get_random_field()``!
+
+The application of any operator is done by one of the :py:meth:`*times` methods that each invoke :py:meth:`_briefing`, :py:meth:`_*multiply` and :py:meth:`_debriefing` consecutively. The briefing and debriefing are generic methods in which in- and output are checked. This check includes a check of the field's properties, weighting of the field values and may include a basis transformation if required. The :py:meth:`_multiply` method, being the concrete part, then carries out the actual computation.
+
+An operator thus interpretes a given input gently, and accepts scalars, lists, arrays (all raised to fields), and fields in the correct conjugate basis (due to the briefing). In turn, the operator tries to restore the input basis in the output field (in the debriefing).
+
+>>> x_space = hp_space(32)
+>>> k_space = x_space.get_codomain()
+>>> f = field(x_space, random="uni", target=k_space)
+>>> I = identity(x_space) # I.domain == x_space
+>>> I(f).domain == x_space
+True
+>>> I(f).target == k_space
+True
+>>> I(f.transform()) # == I.times(f.transform().transform()).transform()
+<nifty.field>
+>>> I(f.transform()).domain == k_space
+True
+>>> I(f.transform()).target == x_space
+True
+>>> I.domain == x_space
+True
+
+.. note:: Transformations between conjugate spaces are performed automatically if required.
+
+Moreover, every operator instance carries three flags, ``uni`` specifying whethter it is unitary, ``sym`` specifying whether it is symmetric, and ``imp`` specifying whethter the correct weighting is already applied internally, otherwise this is done in the (de)briefing.
+
+Diagonal operators
+,,,,,,,,,,,,,,,,,,
+
+Instances of the :py:class:`diagonal_operator` class are initialized
+by providing a domain and a diagonal.
+
+>>> D = diagonal_operator(point_space(8), diag=[2,2,2,2,4,4,4,8])
+>>> print(D)
+<nifty.diagonal_operator>
+>>> D.diag()
+array([ 2.,  2.,  2.,  2.,  4.,  4.,  4.,  8.])
+>>> D.dim()
+array([8, 8])
+
+However, it sometimes comes in handy to promote fields to diagonal operators, vice versa, or to project out the diagonal. Those casts are implemented by :py:meth:`hat` and :py:meth:`hathat`.
+
+>>> D = diagonal_operator(rg_space(10), diag=2)
+>>> D
+<nifty.diagonal_operator>
+>>> D.hat() # (only) diagonal to field
+<nifty.field>
+>>> D.hat().hat() # (only) diagonal to field to diagonal
+<nifty.diagonal_operator>
+>>> D.hathat() # (only) diagonal to diagonal
+<nifty.diagonal_operator>
+>>> D.inverse_hat().val # (only) inverse diagonal to field
+array([ 0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5,  0.5])
+
+The instance method :py:meth:`diag` returns the diagonal at any cost. For non-diagonal operators and operators that are diagonal in another basis this is done by :py:class:`diagonal_probing`, which is not always exact.
+
+.. warning:: The default ``nrun`` for :py:class:`probing` is meant for test runs, but is in most cases insufficient for getting accurate results.
+
+>>> x_space = hp_space(32)
+>>> k_space = x_space.get_codomain()
+>>> D = diagonal_operator(x_space, diag=2) # D.domain == x_space
+>>> D.diag(domain=k_space, nrun=100) # probing in k_space
+array([ 1.99890666 +0.00000000e+00j,  2.00066400 +0.00000000e+00j,
+        1.99810763 +0.00000000e+00j, ...,  2.00000000 -1.87579648e-11j,
+        2.00000000 -1.24384300e-10j,  2.00000000 -2.84119735e-11j])
+
+Assuming that :py:class:`diagonal_operator` is a covariance a zero-mean Gaussian random field can be drawn by :py:meth:`diagonal_operator.get_random_field`. (Confer :ref:`Gaussian random fields <Garfields>`.)
+
+Linear operators
+,,,,,,,,,,,,,,,,
+
+.. currentmodule:: nifty.nifty_explicit
+
+Most of the operator classes in NIFTY describe linear operators that could be represented by an explicit matrix. In practice, an  implicit description is often more efficient, though.
+
+An implicit operator can be cast into an explicit one using the :py:func:`explicify` function that returns an instance of the :py:class:`explicit_operator` class.
+
+>>> space = point_space(3)
+>>> D = diagonal_operator(space, diag=[1, 2, 3])
+>>> E = explicify(D)
+>>> E.get_matrix()
+array([[ 1.,  0.,  0.],
+       [ 0.,  2.,  0.],
+       [ 0.,  0.,  3.]])
+
+The ("bare") matrix is stored as a two-dimensional `numpy.ndarray <http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html>`_ and can be accessed by the ``val`` property, by indexing the field with ``[:,:]``, or  -- as recommended -- by calling :py:meth:`explicit_operator.get_matrix`.
+
+.. note:: The :py:meth:`get_matrix` method takes the keyword argument ``bare``. Compare the "bare" and "non-bare" matrices in case of a continuos space!
+
+.. currentmodule:: nifty
+
+Building your own operator
+,,,,,,,,,,,,,,,,,,,,,,,,,,
+
+Any self-built operator class should be inherited from the :py:class:`operator`. So it inherits the correct substructure and auxiliary methods. Then (almost) all it needs is the development of the :py:meth:`_multiply` method and a proper initialization.
+
+Here, a *non-linear* operator that returns the componentwise squared field is constructed.
+
+>>> class my_operator(operator):
+...
+...     def _multiply(self, x):
+...         return x[:] ** 2
+...
+
+When an instance of ``my_operator`` is created, all relevant properties need to be set. Here, the domain and the ``imp`` flag suffice, and this keeps the operator useable in any space. (Since this is a componentwise operation and does not involve position integrals, the ``imp`` flag is set to ``True``.)
+
+.. warning:: Setting the ``imp`` flag inaccurate can cause faulty algorithmic results that are hard to debug.
+
+>>> space = point_space(8)
+>>> M = my_operator(space, imp=True)
+>>> M(3)
+<nifty.field>
+>>> M(3).val
+array([ 9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.])
+
+>>> space = rg_space(10)
+>>> M = my_operator(space, imp=True)
+>>> M(3).val
+array([ 9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.])
+
+In order to extend the ``my_operator`` class to be of more general use, the ``para`` property can by used.
+
+>>> class my_operator(operator):
+...
+...     def _multiply(self, x):
+...         return x[:] ** self.para
+...
+
+Since there are no constraints on ``para``, it can be of any desired type. (This is especially useful, if several instanced operators need to be combined. List them as one parameter and apply them insinde :py:meth:`_multiply` as required!)
+
+>>> space = rg_space(10)
+>>> M = my_operator(space, imp=True, para=2)
+>>> M(3).val
+array([ 9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.,  9.])
+>>> M = my_operator(space, imp=True, para=np.arange(space.dim()))
+>>> M(2).val
+array([   1.,    2.,    4.,    8.,   16.,   32.,   64.,  128.,  256.,  512.])
+
+`(Read more about operators.) <operator.html>`_
+
+Power & projection operators **(advanced)**
+,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
+
+Although this section is restricted to statistically homogeneous and isotropic Gaussian random fields and power spectra thereof, everything discussed here is generalizable. (Confer :ref:`power indexing <power_index>` and :ref:`Gaussian random fields <Garfields>`.)
+
+The :py:class:`power_operator` is a special :py:class:`diagonal_operator` that handles power spectra properly. A :py:class:`power_operator` performs applications invoked by :py:meth:`*times` analogous to a :py:class:`diagonal_operator`. The difference lies in its initialization, which requires a power spectrum instead of a diagonal, and in the :py:meth:`power_operator.get_random_field` and :py:meth:`power_operator.get_projection_operator` method.
+
+>>> x_space = hp_space(32)
+>>> k_space = x_space.get_codomain()
+>>> spectrum = (lambda k: (1 / (k + 1)) ** 3)
+>>> S = power_operator(k_space, spec=spectrum)
+>>> f = S.get_random_field(domain=x_space) # == field(x_space, target=k_space, random="syn", spec=spectrum)
+>>> f.plot() # a pylab figure will open.
+>>> f.plot(power=True, other=spectrum, mono=False) # another pylab figure will open.
+
+The advantages of a :py:class:`power_operator` are the automated indexing, the methods mentioned above, and the inherited ability to perform basis transformation on fields it is applied to. These concepts become more evident when considering an explicit example.
+
+.. math::
+
+    S\:\big( \phantom{\big|} \vec{k} \:,\: \vec{l} \phantom{\big|} \big) \quad=\quad P\big( \big| \vec{k} \big| \big) \; \delta\big( \vec{k} - \vec{l} \big)
+    \qquad\qquad\qquad
+    S\:\big( \phantom{\big|} \vec{x} \:,\: \vec{y} \phantom{\big|} \big) \quad=\quad S\:\big( \big| \vec{x} - \vec{y} \big| \big)
+
+>>> x_space = rg_space(10, zerocenter=False)
+>>> k_space = x_space.get_codomain()
+>>> S = power_operator(k_space, spec=(lambda k: (1 / (k + 1) ** 2))) # implicit operator
+>>> Skl = explicify(S) # matrix in k_space
+>>> Skl.plot(cmap=pl.cm.gray_r) # a pylab figure will open.
+>>> Sxy = explicify(S, newdomain=x_space) # matrix in x_space
+>>> Sxy.plot(vmin=0, vmax=2, cmap=pl.cm.gray_r) # another pylab figure will open.
+
++---------------------------+---------------------------+
+| .. image:: images/Skl.png | .. image:: images/Sxy.png |
+|     :width:  50 %         |     :width:  50 %         |
++---------------------------+---------------------------+
+
+The :py:class:`projection_operator` class was designed to create
+arbitrary projections. Projections are performed onto a basis subset
+-- the basis is defined by the space specified as domain, and the
+subsets are constructed according to the ``assign`` keyword. (Using
+:py:meth:`power_operator.get_projection_operator` those specifications
+are superfluous, and projections will reproduce the harmonic modes.)
+The initialized :py:meth:`projection_operator` by default projects to
+all specified bands at once (resulting in a nested field), but
+individual projections are possible by giving the keyword argument ``band``.
+
+>>> x_space = hp_space(32)
+>>> k_space = x_space.get_codomain()
+>>> spectrum = (lambda k: (1 / (k + 1)) ** 3)
+>>> S = power_operator(k_space, spec=spectrum)
+>>> f = S.get_random_field(domain=x_space) # == field(x_space, target=k_space, random="syn", spec=spectrum)
+>>> Sk = S.get_projection_operator() # == projection_operator(k_space, assign=k_space.power_indices["pindex"])
+>>> len(k_space.power_indices["kindex"])
+96
+>>> Sk.bands()
+96
+>>> Sk(f)
+<nifty.field>
+>>> print(Sk(f).domain)
+nifty.nested_space instance
+- nest = [<nifty.point_space>, <nifty.hp_space>]
+>>> print(Sk(f).domain.nest[0])
+nifty.point_space instance
+- num      = 96
+- datatype = numpy.float64
+>>> print(Sk(f).domain.nest[1])
+nifty.hp_space instance
+- nside = 32
+>>> Sk(f,band=1)
+<nifty.field>
+>>> print(Sk(f,band=1).domain)
+nifty.hp_space instance
+- nside = 32
+
+How does this look in practice? Below, the quadrupole of the "Faraday Map" is projected out.
+
+>>> x_space = hp_space(128)
+>>> k_space = x_space.get_codomain()
+>>> f = field(x_space, val=np.load("/SOMEWHERE/nifty/demos/demo_faraday_map.npy")) # field values with some easily visible features
+>>> Sk = projection_operator(k_space)
+>>> from nifty.nifty_cmaps import *
+>>> f.plot(vmin=-4, vmax=4, cmap=ncmap.fm()) # a pylab figure will open.
+>>> Sk(f,band=2).plot(vmin=-4, vmax=4, cmap=ncmap.fm()) # another pylab figure will open.
+
++----------------------------+----------------------------+
+| .. image:: images/f_00.png | .. image:: images/f_02.png |
+|     :width:  66 %          |     :width:  66 %          |
++----------------------------+----------------------------+
+
+Do you want to know more?
+.........................
+
+As a continuation of this tutorial, you can browse through the individual class documentations or have a look at the `Demos <demo.html>`_ provided within the NIFTY package.
+
+