Repository: GBillotey/Fractalshades (language: python)

Issue:
Auto potential params
test


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index c299fa8c..97fb42bc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,12 +6,22 @@ dist/
 dev/
 dev_pyQt/
 *.egg-info
-/tests/image_comparison
 *plot_multiproc*
 *multiproc_calc*
 *_local_test*
 *_temp*
+
+# ignore temporary files in tests - keep REFERENCE_DATA !
+/tests/image_comparison
 /tests/known_fails/_data/
 tests/_images_comparison/
-examples/*
-\!examples/*.py
+
+# ignore the files in examples folders except those matching extensions
+examples/**/*.*
+!examples/**/*.py
+!examples/**/*.rst
+
+# ignore the doc build files
+docs/examples/
+docs/_build/
+
diff --git a/dev_notes.txt b/dev_notes.txt
index 9a8e085d..5cf1f25e 100644
--- a/dev_notes.txt
+++ b/dev_notes.txt
@@ -200,6 +200,8 @@ my_struct = types.Record.make_c_struct([
 ])
 
 
+Add an example with :
+http://paulbourke.net/fractals/triangle/
 
 
 
diff --git a/docs/API.rst b/docs/API.rst
new file mode 100644
index 00000000..58b6a643
--- /dev/null
+++ b/docs/API.rst
@@ -0,0 +1,17 @@
+
+Application Programming Interface
+*********************************
+
+Ins this section you will find the detailed documentation of the API:
+
+.. toctree::
+    :maxdepth: 2
+    :caption: Contents:
+
+    API/core.rst
+    API/models.rst
+    API/postproc.rst
+    API/layers.rst
+    API/colors.rst
+
+
diff --git a/docs/API/colors.rst b/docs/API/colors.rst
new file mode 100644
index 00000000..c4cfa735
--- /dev/null
+++ b/docs/API/colors.rst
@@ -0,0 +1,17 @@
+Color mapping
+~~~~~~~~~~~~~
+
+.. module:: fractalshades.colors
+
+This section describe the objects provided to map real areas to colors.
+
+`Fractal_colormap`/s are usually defined through the GUI interface (which
+can export them to python readable format).
+
+For convenience, a `classic_colormap` is provided (blue and orange-brown
+tones), a few examples of the :doc:`../examples/index` use it.
+
+.. autoclass:: Fractal_colormap
+    :members: __init__
+
+.. autodata:: classic_colormap
diff --git a/docs/API/core.rst b/docs/API/core.rst
new file mode 100644
index 00000000..c5fe7c7b
--- /dev/null
+++ b/docs/API/core.rst
@@ -0,0 +1,18 @@
+Core components
+~~~~~~~~~~~~~~~
+
+.. module:: fractalshades
+
+.. autoclass:: Fractal
+    :members: __init__, zoom, run, clean_up
+
+.. autoclass:: PerturbationFractal
+    :members:  __init__, zoom, run
+
+.. autoclass:: Fractal_plotter
+    :members:  __init__, add_layer, plot
+
+.. autodata:: zoom_options
+
+.. autodata:: calc_options
+
diff --git a/docs/API/layers.rst b/docs/API/layers.rst
new file mode 100644
index 00000000..c5f1dc74
--- /dev/null
+++ b/docs/API/layers.rst
@@ -0,0 +1,31 @@
+Image layers
+~~~~~~~~~~~~
+
+Fractalshades graphical export is based on the concept of layers.
+Each layer is linked to a postprocessing field, it can be directly plotted
+or can be used to modify another layer.
+
+.. module:: fractalshades.colors.layers
+
+.. autoclass:: Virtual_layer
+    :members: __init__, set_mask
+
+.. autoclass:: Bool_layer
+    :members:  __init__
+
+.. autoclass:: Color_layer
+    :members:  __init__, set_mask, shade, overlay, set_twin_field 
+
+.. autoclass:: Normal_map_layer
+
+.. autoclass:: Grey_layer
+    :members:  __init__, set_mask
+
+.. autoclass:: Disp_Layer
+    :members:
+
+.. autoclass:: Blinn_lighting
+    :members: __init__, add_light_source
+
+.. autoclass:: Overlay_mode
+    :members: __init__
diff --git a/docs/API/models.rst b/docs/API/models.rst
new file mode 100644
index 00000000..e65c1622
--- /dev/null
+++ b/docs/API/models.rst
@@ -0,0 +1,14 @@
+Fractal models
+~~~~~~~~~~~~~~
+
+This section describes the concrete implementation of
+`fractalshades.Fractal` available.
+
+.. module:: fractalshades.models
+
+.. autoclass:: fractalshades.models.Mandelbrot
+    :members: __init__, zoom, base_calc, newton_calc
+
+.. autoclass:: fractalshades.models.Perturbation_mandelbrot
+    :members:  __init__, zoom, calc_std_div, calc_fast
+
diff --git a/docs/API/postproc.rst b/docs/API/postproc.rst
new file mode 100644
index 00000000..4279995e
--- /dev/null
+++ b/docs/API/postproc.rst
@@ -0,0 +1,40 @@
+Postprocessing fields
+~~~~~~~~~~~~~~~~~~~~~
+
+To make the link between the raw calculation data and the colored layers,
+generic postprocessing are availables. This section describes them.
+
+
+.. module:: fractalshades.postproc
+
+
+.. autoclass:: Postproc
+
+.. autoclass:: Raw_pp
+    :members:  __init__
+
+.. autoclass:: Continuous_iter_pp
+    :members:  __init__
+
+.. autoclass:: DEM_normal_pp
+    :members:  __init__
+
+.. autoclass:: DEM_pp
+    :members:  __init__
+
+.. autoclass:: Attr_normal_pp
+    :members:  __init__
+
+.. autoclass:: Attr_pp
+    :members:  __init__
+
+.. autoclass:: Postproc_batch
+    :members: __init__, add_postproc
+
+.. autoclass:: Fractal_array
+    :members:  __init__
+
+
+
+
+
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 00000000..298ea9e2
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,19 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+SOURCEDIR     = .
+BUILDDIR      = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+	@$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option.  $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+	@$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/docs/_static/logo2.jpg b/docs/_static/logo2.jpg
new file mode 100644
index 00000000..44d9d871
Binary files /dev/null and b/docs/_static/logo2.jpg differ
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 00000000..0f0cd806
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,242 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# 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
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+import os
+import sys
+import fractalshades
+fs_dir = os.path.abspath(os.path.dirname(fractalshades.__file__))
+sys.path.insert(0, fs_dir)
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'Fractalshades'
+copyright = '2021, Geoffroy Billotey'
+author = 'Geoffroy Billotey'
+
+# The short X.Y version
+version = fractalshades.__version__
+# The full version, including alpha/beta/rc tags
+release = fractalshades.__version__
+
+
+# -- General configuration ---------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+    'sphinx.ext.autodoc',
+    'sphinx.ext.autosummary',
+    # 'sphinx.ext.doctest',
+    # 'sphinx.ext.inheritance_diagram',
+    'sphinx.ext.imgmath',
+    'sphinx.ext.viewcode',
+    # 'sphinx.ext.napoleon',  # for NumPy style docstrings
+    'sphinx.ext.githubpages',
+    'numpydoc',  # Needs to be loaded *after* autodoc.
+    'sphinx_gallery.gen_gallery',
+    #'sphinxext.github',
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx' #None
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'sphinx_rtd_theme' # 'alabaster'
+
+# 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 static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# The default sidebars (for documents that don't match any pattern) are
+# defined by theme itself.  Builtin themes are using these templates by
+# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
+# 'searchbox.html']``.
+#
+# html_sidebars = {}
+
+
+# -- Options for HTMLHelp output ---------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'Fractalshadesdoc'
+
+
+# -- Options for LaTeX output ------------------------------------------------
+
+latex_elements = {
+    # The paper size ('letterpaper' or 'a4paper').
+    #
+    # 'papersize': 'letterpaper',
+
+    # The font size ('10pt', '11pt' or '12pt').
+    #
+    # 'pointsize': '10pt',
+
+    # Additional stuff for the LaTeX preamble.
+    #
+    # 'preamble': '',
+
+    # Latex figure (float) alignment
+    #
+    # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+#  author, documentclass [howto, manual, or own class]).
+latex_documents = [
+    (master_doc, 'Fractalshades.tex', 'Fractalshades Documentation',
+     'Geoffroy Billotey', 'manual'),
+]
+
+
+# -- Options for manual page output ------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+    (master_doc, 'fractalshades', 'Fractalshades Documentation',
+     [author], 1)
+]
+
+
+# -- Options for Texinfo output ----------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+#  dir menu entry, description, category)
+texinfo_documents = [
+    (master_doc, 'Fractalshades', 'Fractalshades Documentation',
+     author, 'Fractalshades', 'One line description of project.',
+     'Miscellaneous'),
+]
+
+
+# -- Options for Epub output -------------------------------------------------
+
+# Bibliographic Dublin Core info.
+epub_title = project
+
+# The unique identifier of the text. This can be a ISBN number
+# or the project homepage.
+#
+# epub_identifier = ''
+
+# A unique identification for the text.
+#
+# epub_uid = ''
+
+# A list of files that should not be packed into the epub file.
+epub_exclude_files = ['search.html']
+
+
+# -- Extension configuration -------------------------------------------------
+
+# The base URL which points to the root of the HTML documentation.
+# It is used to indicate the location of document using The Canonical Link
+# Relation. Default: ''.
+# Used by also by sphinx.ext.githubpages CNAME
+html_baseurl = ""
+html_logo = "_static/logo2.jpg"
+
+default_role = 'py:obj'
+
+# Not sure how this works... to be investigated
+# extlinks = {'numpy': ('https://numpy.org/', '%s')}
+
+# numpydoc config
+numpydoc_show_class_members = False
+
+# Automatic build of the gallery
+import sphinx_gallery
+# We will not write to disk but feed fractalshades_scraper
+fractalshades.settings.output_mode = "Pillow" 
+
+# https://sphinx-gallery.github.io/stable/advanced.html?highlight=scrapers#example-2-detecting-image-files-on-disk
+def fractalshades_scraper(block, block_vars, gallery_conf):
+    # We use a list to collect references to image names
+    image_names = list()
+    # The `image_path_iterator` is created by Sphinx-Gallery, it will yield
+    # a path to a file name that adheres to Sphinx-Gallery naming convention.
+    image_path_iterator = block_vars['image_path_iterator']
+
+    # Define a list of our already-created figure objects.
+    list_of_my_figures = fractalshades.get_figures()
+
+    # Iterate through figure objects, save to disk, and keep track of paths.
+    for fig, image_path in zip(list_of_my_figures, image_path_iterator):
+        fig.save_png(image_path)
+        image_names.append(image_path)
+
+    # Close all references to figures so they aren't used later.
+    fractalshades.close('all')
+
+    # Use the `figure_rst` helper function to generate the rST for this
+    # code block's figures. Alternatively you can define your own rST.
+    return sphinx_gallery.scrapers.figure_rst(image_names,
+                                              gallery_conf['src_dir'])
+
+sphinx_gallery_conf = {
+     'examples_dirs': '../examples',   # path to your example scripts
+     'gallery_dirs': 'examples',       # path to where to save gallery generated output
+     'filename_pattern': r'\.py',     # all python files r'\.py',
+     'ignore_pattern': r'sgskip\.py', # except those ending with sgskip
+     'image_scrapers': (fractalshades_scraper),
+     # 'plot_gallery': 'False',  # Activate this to skip the calculation
+}
diff --git a/docs/extending.rst b/docs/extending.rst
new file mode 100644
index 00000000..80409f04
--- /dev/null
+++ b/docs/extending.rst
@@ -0,0 +1,110 @@
+Extending
+*********
+
+This section provides details for the advanced used on how to 
+extend fractalshades.
+
+new coloring technique
+~~~~~~~~~~~~~~~~~~~~~~
+
+To add a new coloring technique, one should first consider subclassing
+`fractalshades.postproc.Postproc`. If the needed raw complex or integer
+fields are already available, this will be enough.
+
+If a new calculation output is needed, it will be necessary to also add a
+new calculation routine to one of the `fractalshades.Fractal` models.
+This would be the case for instance if you wish to add a custom orbit trap.
+This can either be done through sub-classing or monkey-patching. Remember to
+tag (decorate) your routine with @ `fractalshades.calc_options`.
+
+new fractal
+~~~~~~~~~~~
+
+To define a new fracal, subclass `fractalshades.Fractal`. This is almost
+straightforward for standard (double-precision) fractals as one only has to 
+follow the steps :
+
+- Defines the list of complex codes (corresponding to the complex fields
+  passed for each point at each iteration), integer codes (the integer fields)
+  and stop codes (the reason for exiting the loop) ;
+  store them as class attributes
+
+  .. code-block:: python
+
+      self.codes = (complex_codes, int_codes, stop_codes)
+
+
+- Initialise the datatypes, usually with call to `init_data_types`:
+
+  .. code-block:: python
+
+      self.init_data_types(np.complex128)
+
+- Define the initialisation function and store it as a class attribute (Note
+  that we use a wrapper to allow customization). This step is done with 
+  classic numpy calls.
+  As an example if we only have 4 complex fields
+  [“zn”, “dzndz”, “dzndc”, “d2zndc2”] and no integer field:
+
+  .. code-block:: python
+
+        def initialize():
+            def func(Z, U, c, chunk_slice):
+                Z[3, :] = 0.
+                Z[2, :] = 0.
+                Z[1, :] = 0.
+                Z[0, :] = 0.
+            return func
+        self.initialize = initialize
+
+- Define the iteration function and store it as a class attribute (Note
+  that we use a wrapper to allow customization). This step is time
+  critical so is done with numba jitted functions (njit).
+
+  As an example if we only have 4 complex fields 
+  ["zn", "dzndz", "dzndc", "d2zndc2"] and no integer field:
+
+  .. code-block:: python
+
+        def iterate():
+            Mdiv_sq = self.M_divergence ** 2
+            epscv_sq = self.epsilon_stationnary ** 2
+            max_iter = self.max_iter
+            @numba.njit
+            def numba_impl(Z, U, c, stop_reason, n_iter):
+                if n_iter >= max_iter:
+                    stop_reason[0] = 0
+
+                Z[3] = 2 * (Z[3] * Z[0] + Z[2]**2)
+                Z[2] = 2 * Z[2] * Z[0] + 1.
+                Z[1] = 2 * Z[1] * Z[0] 
+                Z[0] = Z[0]**2 + c
+                if n_iter == 1:
+                    Z[1] = 1.
+
+                if Z[0].real ** 2 + Z[0].imag ** 2 > Mdiv_sq:
+                    stop_reason[0] = 1
+                if Z[1].real ** 2 + Z[1].imag ** 2 < epscv_sq:
+                    stop_reason[0] = 2
+
+            return numba_impl
+        self.iterate = iterate
+
+- Remember that inner loop will exit whenever stop_reason[0] is set to a
+  value != -1. The value shall be a null or positive integer : 
+  corresponding to the index in `stop_codes`. For instance with `stop_codes`
+  ["max_iter", "divergence", "stationnary"] if you reach the criteria for
+  divergence you have noticed line :
+
+  .. code-block:: python
+
+        if Z[0].real ** 2 + Z[0].imag ** 2 > Mdiv_sq:
+            stop_reason[0] = 1
+
+Et voila (seee also `fractalshades.models.Mandelbrot.base_calc` source code).
+To be adapted for your needs...
+
+For fractals implementing perturbation technique, this is beyong the scope
+of this use manual.
+
+
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
new file mode 100644
index 00000000..17c6dc89
--- /dev/null
+++ b/docs/getting_started.rst
@@ -0,0 +1,85 @@
+
+Getting started
+***************
+
+Installing
+~~~~~~~~~~
+
+Fractalshades is published on the Python Package Index (PyPI_). To install the
+last published version and its dependencies, you should run: [#f1]_
+
+.. code-block:: console
+
+    python3 -m pip install --user fractalshades
+
+or, to install directly the latest version from Github master:
+
+.. code-block:: console
+
+    python -m pip install https://github.com/GBillotey/Fractalshades.git
+
+This package relies on the following dependencies:
+
+- numpy_
+- numba_
+- mpmath_
+- gmpy2_
+- Pillow_
+- PyQt5_
+
+.. _numpy: https://numpy.org/
+.. _numba: http://numba.pydata.org/
+.. _mpmath: https://mpmath.org/
+.. _gmpy2: https://gmpy2.readthedocs.io/en/latest/
+.. _Pillow: https://pillow.readthedocs.io/en/stable/
+.. _PyQt5: https://pypi.org/project/PyQt5/
+.. _PyPI: https://pypi.org/
+
+They should install automatically through `pip`. A special case is gmpy2 as it
+needs the most recent versions of GMP, MPFR and MPC multi-precision
+arithmetic libraries. If your distribution does not include them you will have
+to install them manually. For instance under Linux:
+
+.. code-block:: console
+
+    sudo apt-get install libgmp-dev
+    sudo apt-get install libmpfr-dev
+    sudo apt-get install libmpc-dev
+
+.. [#f1] These instructions describe installation to your Python home
+         directory. You could also consider the installation of
+         `Fractalshades` in a virtual environment (a self-contained directory
+         tree that contains a Python installation for a particular version of
+         Python), through venv_.
+
+.. _venv: https://docs.python.org/3/tutorial/venv.html
+
+
+
+A 5-minutes guide to fractalshades
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Fractalshades exposes 3 kinds of component each implementing a different
+functionality:
+
+  - The core or calculation components runs under the hood all the
+    calculations necessary for a plot and
+    stores the intermediate raw results. This is typically done in a subclass
+    of `fractalshades.Fractal` base class.
+    The list of the available fractal models is found here :
+    :doc:`API/models`.
+
+  - The plotting components will open the raw results and apply user-selected
+    post-processing, to generate the image output. The base class for this
+    part is `fractalshades.Fractal_plotter`.
+    Common post-processing routines are available, they are listed under
+    the :doc:`API/postproc` section.
+
+  - In order to explore a fractal and select a location, a GUI is necessary.
+    Fractalshades comes with a very flexible and user-configurable graphical
+    interface based on PyQt5 framework.
+
+The best way to start is probably to have a look at the 
+:doc:`examples/index` section.
+
+
diff --git a/docs/index.rst b/docs/index.rst
new file mode 100644
index 00000000..760ad13e
--- /dev/null
+++ b/docs/index.rst
@@ -0,0 +1,31 @@
+
+.. title:: Fractalshades documentation
+
+.. module:: fractalshades
+
+
+Fractalshades |version| documentation
+=====================================
+
+Fractalshades is a library for creating static and interactive visualizations 
+or Fractals in Python.
+
+
+.. toctree::
+    :maxdepth: 2
+    :caption: Contents:
+
+    overview.rst
+    getting_started.rst
+    examples/index.rst
+    API.rst
+    extending.rst
+
+
+Indices and tables
+******************
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
+
diff --git a/docs/overview.rst b/docs/overview.rst
new file mode 100644
index 00000000..349abf38
--- /dev/null
+++ b/docs/overview.rst
@@ -0,0 +1,27 @@
+Overview
+********
+
+Fractalshades is a library for creating static and interactive visualizations 
+or fractals in Python. It targets primarily the escape-time fractals, a wide 
+class of fractals among which the most famous is probably the Mandelbrot_ set.
+
+.. _Mandelbrot: https://en.wikipedia.org/wiki/Mandelbrot_set
+
+Fractalshades aims at:
+
+  - providing built-in formulas for most common escape-time fractal, which 
+    let you produce beautiful images with just a few lines of code
+
+  - for more advanced users, providing easy ways of implementing custom 
+    formulas or coloring options, which will then benefit of the same 
+    time-optimisation and graphical user interface.
+
+Some fractals are only available at double precision, while some also allow
+ultra-deep zooming (1e-300 scale and beyong) thanks to state-of-the-art
+algorithm (based on perturbation technique).
+
+The time-critical core-loops benefit from multiprocesssing and leverage
+just-in-time compiling through numba.
+
+
+
diff --git a/docs/readme.txt b/docs/readme.txt
new file mode 100644
index 00000000..547cbc53
--- /dev/null
+++ b/docs/readme.txt
@@ -0,0 +1,13 @@
+To build the docs you will need to install sphinx, sphinx_rtd_theme etc...
+
+
+> python3 -m pip install --user sphinx
+> python3 -m pip install --user sphinx_rtd_theme
+> python3 -m pip install --user sphinx-gallery
+
+
+Then (under Linux) from a command prompt in this directory :
+> make html 
+
+
+https://sphinx-gallery.github.io/stable/getting_started.html
diff --git a/examples/README.rst b/examples/README.rst
new file mode 100644
index 00000000..693f4ac6
--- /dev/null
+++ b/examples/README.rst
@@ -0,0 +1,12 @@
+*******
+Gallery
+*******
+
+This gallery presents a showcase of the possibilities of fractalshades.
+The images are auto-generated from the exemples source code, to keep both the
+running time and the size of the documentation acceptable we limit the images
+quality here. Download and run the file with higher quality at home !
+
+Click on an image to see the source code
+
+
diff --git a/examples/batch_mode/README.rst b/examples/batch_mode/README.rst
new file mode 100644
index 00000000..22436e4d
--- /dev/null
+++ b/examples/batch_mode/README.rst
@@ -0,0 +1,3 @@
+Batch mode gallery
+==================
+
diff --git a/examples/batch_mode/double_embedded_julia.py b/examples/batch_mode/double_embedded_julia.py
new file mode 100644
index 00000000..9bb6c19c
--- /dev/null
+++ b/examples/batch_mode/double_embedded_julia.py
@@ -0,0 +1,154 @@
+# -*- coding: utf-8 -*-
+"""
+=========================
+Double embedded Julia set
+=========================
+Example plot a double-embedded Julia set in mandelbrot power-2
+
+This example even if not too deep is beyong the reach of double precision
+[#f1]_, at *7e-22*. We will use a `Perturbation_mandelbrot` instance.
+
+
+.. [#f1] **Credit:** Coordinates from Robert P. Munafo website:
+        <https://mrob.com/pub/muency/secondorderembeddedjuliase.html>
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    DEM_normal_pp,
+    Fieldlines_pp,
+    Raw_pp,
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+    Normal_map_layer,
+    Virtual_layer,
+    Blinn_lighting,
+)
+
+
+def plot(plot_dir):
+    # A simple showcase using perturbation technique
+    x = "-1.768667862837488812627419470"
+    y = "0.001645580546820209430325900"
+    dx = "7.0e-22"
+    precision = 30
+    nx = 800
+
+    calc_name="mandelbrot"
+    
+    c3 = np.array([125, 225, 255]) / 255.
+    c0 = np.array([1, 254, 255]) / 255.
+    c1 = c0 * 0.01
+
+    colors = np.vstack((c0[np.newaxis, :],
+                        c1[np.newaxis, :],
+                        c3[np.newaxis, :]))
+    colormap = fscolors.Fractal_colormap(kinds="Lch", colors=colors,
+         grad_npts=200, grad_funcs="x**0.5", extent="mirror")
+
+    # Run the calculation
+    f = fsm.Perturbation_mandelbrot(plot_dir)
+    f.zoom( precision=precision,
+            x=x,
+            y=y,
+            dx=dx,
+            nx=nx,
+            xy_ratio=1.0,
+            theta_deg=0., 
+            projection="cartesian",
+            antialiasing=False)
+
+    f.calc_std_div(
+            datatype=np.complex128,
+            calc_name=calc_name,
+            subset=None,
+            max_iter=200000,
+            M_divergence=1.e3,
+            epsilon_stationnary=1.e-3,
+            SA_params={"cutdeg": 64,
+                       "cutdeg_glitch": 8,
+                       "SA_err": 1.e-4},
+            glitch_eps=1.e-6,
+            interior_detect=False,
+            glitch_max_attempt=2)
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    pp.add_postproc("DEM_map", DEM_normal_pp(kind="potential"))
+    pp.add_postproc("fieldlines",
+                Fieldlines_pp(n_iter=8, swirl=0.4, damping_ratio=0.1))
+
+    plotter = fs.Fractal_plotter(pp)   
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Normal_map_layer("DEM_map", max_slope=60, output=False))
+    plotter.add_layer(Virtual_layer("fieldlines", func=None, output=False))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            #probes_z=[0.095, 0.11],
+            probes_z=[0.07, 0.075],
+            probes_kind="absolute",
+            output=True
+    ))
+
+    plotter["cont_iter"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+    plotter["DEM_map"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+
+    # This is the line where we indicate that coloring is a combination of
+    # "Continuous iteration" and "fieldines values"
+    plotter["cont_iter"].set_twin_field(plotter["fieldlines"], 0.0001)
+
+    plotter.plot()
+    # This is where we define the lighting (here 3 ccolored light sources)
+    # and apply the shading
+    light = Blinn_lighting(0.2, np.array([1., 1., 1.]))
+    light.add_light_source(
+        k_diffuse=0.2,
+        k_specular=200.,
+        shininess=1400.,
+        angles=(45., 20.),
+        coords=None,
+        color=np.array([1.0, 1.0, 0.5]))
+    light.add_light_source(
+        k_diffuse=1.,
+        k_specular=0.,
+        shininess=0.,
+        angles=(-905., 20.),
+        coords=None,
+        color=np.array([1., 1., 0.8]))
+    light.add_light_source(
+        k_diffuse=1.,
+        k_specular=0.,
+        shininess=0.,
+        angles=(-15., 20.),
+        coords=None,
+        color=np.array([1., 1., 0.5]))
+    plotter["cont_iter"].shade(plotter["DEM_map"], light)
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/run_perturb3.py b/examples/batch_mode/run_perturb3_sgskip.py
similarity index 100%
rename from examples/run_perturb3.py
rename to examples/batch_mode/run_perturb3_sgskip.py
diff --git a/examples/run_perturb.py b/examples/batch_mode/run_perturb_sgskip.py
similarity index 100%
rename from examples/run_perturb.py
rename to examples/batch_mode/run_perturb_sgskip.py
diff --git a/examples/batch_mode/seahorse_basic.py b/examples/batch_mode/seahorse_basic.py
new file mode 100644
index 00000000..6cc4ea05
--- /dev/null
+++ b/examples/batch_mode/seahorse_basic.py
@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+"""
+======================
+Seahorse basic example
+======================
+
+This basic example shows how to create a color layer, displaying the 
+"continuous iteration number" for Mandelbrot (power 2) fractal.
+
+The location is a shallow one in the main Seahorse valley.
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    Raw_pp,
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+)
+
+def plot(plot_dir):
+    """
+    A very simple example: full view of the Mandelbrot set with escape-time
+    coloring
+    """
+    # Define the parameters for this calculation
+    x = -0.746223962861
+    y = -0.0959468433527
+    dx = 0.00745
+    nx = 800
+
+    calc_name="mandelbrot"
+    colormap = fscolors.classic_colormap
+
+    # Run the calculation
+    f = fsm.Mandelbrot(plot_dir)
+    f.zoom(x=x, y=y, dx=dx, nx=nx, xy_ratio=1.0,
+           theta_deg=0., projection="cartesian", antialiasing=False)
+    f.base_calc(
+        calc_name=calc_name,
+        subset=None,
+        max_iter=1000,
+        M_divergence=100.,
+        epsilon_stationnary= 0.001,
+        datatype=np.complex128)
+    # f.clean_up(calc_name) 
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    plotter = fs.Fractal_plotter(pp)   
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=True
+    ))
+    plotter["cont_iter"].set_mask(
+            plotter["interior"],
+            mask_color=(0., 0., 0.)
+    )
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/batch_mode/seahorse_fieldlines1.py b/examples/batch_mode/seahorse_fieldlines1.py
new file mode 100644
index 00000000..9149261f
--- /dev/null
+++ b/examples/batch_mode/seahorse_fieldlines1.py
@@ -0,0 +1,97 @@
+# -*- coding: utf-8 -*-
+"""
+===========================================
+Seahorse fieldlines example "tint_or_shade"
+===========================================
+
+This example shows one of the ways to plot fieldlines: here the fieldlines 
+values are used to modify the original layer colors : they are tinted or
+shaded.
+
+The location is a shallow one in the main Seahorse valley.
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    Fieldlines_pp,
+    Raw_pp,
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+    Grey_layer,
+    Overlay_mode
+)
+
+def plot(plot_dir):
+    """
+    Using field lines : a shallow zoom in the Seahorses valley
+    Coloring based on continuous iteration
+    Colors 'tinted or shaded' by the fieldliles
+    """
+    # Define the parameters for this calculation
+    x = -0.746223962861
+    y = -0.0959468433527
+    dx = 0.00745
+    nx = 800
+    calc_name="mandelbrot"
+    colormap = fscolors.classic_colormap
+
+    # Run the calculation
+    f = fsm.Mandelbrot(plot_dir)
+    f.zoom(x=x, y=y, dx=dx, nx=nx, xy_ratio=1.0,
+           theta_deg=0., projection="cartesian", antialiasing=False)
+    f.base_calc(
+        calc_name=calc_name,
+        subset=None,
+        max_iter=5000,
+        M_divergence=100.,
+        epsilon_stationnary= 0.001,
+        datatype=np.complex128)
+    # f.clean_up(calc_name) # keep this line if you want to force recalculate
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    pp.add_postproc("fieldlines",
+                Fieldlines_pp(n_iter=4, swirl=1., damping_ratio=0.1))
+
+    plotter = fs.Fractal_plotter(pp)   
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=True
+    ))
+    plotter.add_layer(Grey_layer("fieldlines", func=None, output=False))
+    plotter["cont_iter"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+    # This is the lines where we indicate that coloring is shaded or tinted
+    # depending on "fieldines" values
+    overlay_mode = Overlay_mode("tint_or_shade", pegtop=1.)
+    plotter["cont_iter"].overlay(plotter["fieldlines"], overlay_mode)
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/batch_mode/seahorse_fieldlines2.py b/examples/batch_mode/seahorse_fieldlines2.py
new file mode 100644
index 00000000..3fdd28c6
--- /dev/null
+++ b/examples/batch_mode/seahorse_fieldlines2.py
@@ -0,0 +1,95 @@
+# -*- coding: utf-8 -*-
+"""
+=======================================
+Seahorse fieldlines example "twinfield"
+=======================================
+
+This example shows one of the ways to plot fieldlines: here the fieldlines 
+values are used to modify the original layer values before applying the
+colormap.
+
+The location is a shallow one in the main Seahorse valley.
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    Fieldlines_pp,
+    Raw_pp,
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+    Virtual_layer
+)
+
+def plot(plot_dir):
+    """
+    Using field lines : a shallow zoom in the Seahorses valley
+    Coloring based on continuous iteration + fieldlines
+    """
+    # Define the parameters for this calculation
+    x = -0.746223962861
+    y = -0.0959468433527
+    dx = 0.00745
+    nx = 800
+    calc_name="mandelbrot"
+    colormap = fscolors.classic_colormap
+
+    # Run the calculation
+    f = fsm.Mandelbrot(plot_dir)
+    f.zoom(x=x, y=y, dx=dx, nx=nx, xy_ratio=1.0,
+           theta_deg=0., projection="cartesian", antialiasing=False)
+    f.base_calc(
+        calc_name=calc_name,
+        subset=None,
+        max_iter=5000,
+        M_divergence=100.,
+        epsilon_stationnary= 0.001,
+        datatype=np.complex128)
+    # f.clean_up(calc_name)
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    pp.add_postproc("fieldlines",
+                Fieldlines_pp(n_iter=4, swirl=1., damping_ratio=0.1))
+
+    plotter = fs.Fractal_plotter(pp)   
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=True
+    ))
+    plotter.add_layer(Virtual_layer("fieldlines", func=None, output=False))
+    plotter["cont_iter"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+
+    # This is the line where we indicate that coloring is a combination of
+    # "Continuous iteration" and "fieldines values"
+    plotter["cont_iter"].set_twin_field(plotter["fieldlines"], 0.15)
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/batch_mode/seahorse_interior.py b/examples/batch_mode/seahorse_interior.py
new file mode 100644
index 00000000..d1c4917b
--- /dev/null
+++ b/examples/batch_mode/seahorse_interior.py
@@ -0,0 +1,177 @@
+# -*- coding: utf-8 -*-
+"""
+=========================
+Seahorse interior example
+=========================
+
+This example shows how to use several layers to plot both the divergent and
+convergent part of a fractal, with fieldlines and scene lighting.
+
+The location is a shallow one in the main Seahorse valley.
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    DEM_normal_pp,
+    Fieldlines_pp,
+    Raw_pp,
+    Fractal_array,
+    Attr_pp,
+    Attr_normal_pp
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+    Normal_map_layer,
+    Grey_layer,
+    Blinn_lighting,
+    Overlay_mode
+)
+
+
+def plot(plot_dir=None):
+    """
+    Combining all : a shallow zoom in the Seahorses valley
+    Coloring of escaping points based on continuous iteration + lighting with a
+    normal maps from distance estimation method
+    Coloring of interior points based on the attracting cycle attractivity
+    """
+    # Define the parameters for this calculation
+    x = -0.746223962861
+    y = -0.0959468433527
+    dx = 0.00745
+    nx = 800
+
+    calc_name="escaping"
+    colormap = fscolors.classic_colormap
+    colormap_int = fscolors.classic_colormap
+
+    # Run the calculation
+    f = fsm.Mandelbrot(plot_dir)
+    f.zoom(x=x, y=y, dx=dx, nx=nx, xy_ratio=1.0,
+           theta_deg=0., projection="cartesian", antialiasing=False)
+    f.base_calc(
+        calc_name=calc_name,
+        subset=None,
+        max_iter=5000,
+        M_divergence=100.,
+        epsilon_stationnary= 0.001,
+        datatype=np.complex128)
+    # f.clean_up("escaping") # keep this line if you want to force recalculate
+    f.run()
+    
+    # Run the calculation for the interior points
+    interior = Fractal_array(f, "escaping", "stop_reason",
+                         func= lambda x: (x != 1))
+    f.newton_calc(
+        calc_name="interior",
+        subset=interior,
+        known_orders=None,
+        max_order=1500,
+        max_newton=20,
+        eps_newton_cv=1.e-12,
+        datatype=np.complex128)
+    # f.clean_up("interior")
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    pp.add_postproc("div", Raw_pp("stop_reason", func="x == 1."))
+    pp.add_postproc("DEM_map", DEM_normal_pp(kind="potential"))
+    pp.add_postproc("fieldlines",
+                Fieldlines_pp(n_iter=10, swirl=1., damping_ratio=0.1))
+    
+    # Defines a second pastproc batch for interior points
+    pp_int = Postproc_batch(f, "interior")
+    pp_int.add_postproc("attr_map", Attr_normal_pp())
+    pp_int.add_postproc("attr", Attr_pp())
+
+    plotter = fs.Fractal_plotter([pp, pp_int])  
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Bool_layer("div", output=False))
+    plotter.add_layer(Normal_map_layer("DEM_map", max_slope=60, output=False))
+    plotter.add_layer(Normal_map_layer("attr_map", max_slope=90, output=False))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=True
+    ))
+    plotter.add_layer(Color_layer(
+            "attr",
+            func=None, #"x * -1. + 1.", #"np.log(x)",
+            colormap=colormap_int,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=False))
+    plotter.add_layer(Grey_layer("fieldlines", func=None,
+                                 curve=lambda x: 0.5 + (x - 0.5) * 0.2,
+                                 output=False))
+
+    # plotter["cont_iter"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+    plotter["DEM_map"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+    plotter["attr"].set_mask(plotter["div"], mask_color=(0.9568, 0.8039, 0.9372))
+
+    # This is where we define the lighting (here 3 ccolored light sources)
+    # and apply the shading
+    light = Blinn_lighting(0.2, np.array([1., 1., 1.]))
+    light.add_light_source(
+        k_diffuse=0.2,
+        k_specular=10.,
+        shininess=400.,
+        angles=(-135., 20.),
+        coords=None,
+        color=np.array([0.05, 0.05, 1.0]))
+    light.add_light_source(
+        k_diffuse=0.2,
+        k_specular=10.,
+        shininess=400.,
+        angles=(135., 20.),
+        coords=None,
+        color=np.array([0.5, 0.5, .4]))
+    light.add_light_source(
+        k_diffuse=1.3,
+        k_specular=0.,
+        shininess=0.,
+        angles=(90., 40.),
+        coords=None,
+        color=np.array([1.0, 1.0, 1.0]))
+
+    # Adds some shading based on the previouly defined normal maps
+    plotter["cont_iter"].shade(plotter["DEM_map"], light)
+    plotter["attr"].shade(plotter["attr_map"], light)
+
+    # Overlay : tint or shade depending on fieldlines layer value
+    overlay_mode = Overlay_mode("tint_or_shade", pegtop=1.)
+    plotter["cont_iter"].overlay(plotter["fieldlines"], overlay_mode)
+
+    # Overlay : alpha composite with "interior" layer ie, where it is not
+    # masked, we take the value of the "attr" layer
+    overlay_mode = Overlay_mode("alpha_composite")
+    plotter["cont_iter"].overlay(plotter["attr"], overlay_mode=overlay_mode)
+
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/batch_mode/seahorse_shaded.py b/examples/batch_mode/seahorse_shaded.py
new file mode 100644
index 00000000..e38c7b16
--- /dev/null
+++ b/examples/batch_mode/seahorse_shaded.py
@@ -0,0 +1,120 @@
+# -*- coding: utf-8 -*-
+"""
+=======================
+Seahorse shaded example
+=======================
+
+This example shows how to create a normal map layer, and link it to a base
+color layer to enable scene lighting.
+
+The location is a shallow one in the main Seahorse valley.
+"""
+
+import os
+import numpy as np
+
+import fractalshades as fs
+import fractalshades.models as fsm
+import fractalshades.colors as fscolors
+from fractalshades.postproc import (
+    Postproc_batch,
+    Continuous_iter_pp,
+    DEM_normal_pp,
+    Raw_pp,
+)
+from fractalshades.colors.layers import (
+    Color_layer,
+    Bool_layer,
+    Normal_map_layer,
+    Blinn_lighting
+)
+
+def plot(plot_dir):
+    """
+    Using lighting : a shallow zoom in the Seahorses valley
+    Coloring based on continuous iteration + lighting with a normal maps from
+    distance estimation method
+    """
+    # Define the parameters for this calculation
+    x = -0.746223962861
+    y = -0.0959468433527
+    dx = 0.00745
+    nx = 800
+
+    calc_name="mandelbrot"
+    colormap = fscolors.classic_colormap
+
+    # Run the calculation
+    f = fsm.Mandelbrot(plot_dir)
+    f.zoom(x=x, y=y, dx=dx, nx=nx, xy_ratio=1.0,
+           theta_deg=0., projection="cartesian", antialiasing=False)
+    f.base_calc(
+        calc_name=calc_name,
+        subset=None,
+        max_iter=5000,
+        M_divergence=100.,
+        epsilon_stationnary= 0.001,
+        datatype=np.complex128)
+    # f.clean_up(calc_name) # keep this line if you want to force recalculation
+    f.run()
+
+    # Plot the image
+    pp = Postproc_batch(f, calc_name)
+    pp.add_postproc("cont_iter", Continuous_iter_pp())
+    pp.add_postproc("interior", Raw_pp("stop_reason", func="x != 1."))
+    pp.add_postproc("DEM_map", DEM_normal_pp(kind="potential"))
+
+    plotter = fs.Fractal_plotter(pp)   
+    plotter.add_layer(Bool_layer("interior", output=False))
+    plotter.add_layer(Normal_map_layer("DEM_map", max_slope=60, output=True))
+    plotter.add_layer(Color_layer(
+            "cont_iter",
+            func="np.log(x)",
+            colormap=colormap,
+            probes_z=[1., 2.],
+            probes_kind="absolute",
+            output=True
+    ))
+
+    plotter["cont_iter"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+    plotter["DEM_map"].set_mask(plotter["interior"], mask_color=(0., 0., 0.))
+
+    # This is where we define the lighting (here 3 ccolored light sources)
+    # and apply the shading
+    light = Blinn_lighting(0.2, np.array([1., 1., 1.]))
+    light.add_light_source(
+        k_diffuse=0.2,
+        k_specular=10.,
+        shininess=400.,
+        angles=(-135., 20.),
+        coords=None,
+        color=np.array([0.05, 0.05, 1.0]))
+    light.add_light_source(
+        k_diffuse=0.2,
+        k_specular=10.,
+        shininess=400.,
+        angles=(135., 20.),
+        coords=None,
+        color=np.array([0.5, 0.5, .4]))
+    light.add_light_source(
+        k_diffuse=1.3,
+        k_specular=0.,
+        shininess=0.,
+        angles=(90., 40.),
+        coords=None,
+        color=np.array([1.0, 1.0, 1.0]))
+    plotter["cont_iter"].shade(plotter["DEM_map"], light)
+    plotter.plot()
+
+
+if __name__ == "__main__":
+    # Some magic to get the directory for plotting: with a name that matches
+    # the file or a temporary dir if we are building the documentation
+    try:
+        realpath = os.path.realpath(__file__)
+        plot_dir = os.path.splitext(realpath)[0]
+        plot(plot_dir)
+    except NameError:
+        import tempfile
+        with tempfile.TemporaryDirectory() as plot_dir:
+            plot(plot_dir)
diff --git a/examples/interactive_mode/README.rst b/examples/interactive_mode/README.rst
new file mode 100644
index 00000000..af32c53c
--- /dev/null
+++ b/examples/interactive_mode/README.rst
@@ -0,0 +1,4 @@
+GUI-mode examples
+=================
+
+Below is a gallery of run with pyQt5
diff --git a/examples/run_interactive.py b/examples/interactive_mode/run_interactives_sgskip.py
similarity index 98%
rename from examples/run_interactive.py
rename to examples/interactive_mode/run_interactives_sgskip.py
index 6855002c..a4041aae 100644
--- a/examples/run_interactive.py
+++ b/examples/interactive_mode/run_interactives_sgskip.py
@@ -122,9 +122,9 @@ def func(fractal: fsm.Perturbation_mandelbrot=fractal,
         fractal.zoom(precision=dps, x=x, y=y, dx=dx, nx=nx, xy_ratio=xy_ratio,
              theta_deg=0., projection="cartesian", antialiasing=False)
 
-        fractal.calc_std_div(complex_type=np.complex128, calc_name=calc_name,
+        fractal.calc_std_div(datatype=np.complex128, calc_name=calc_name,
             subset=None, max_iter=max_iter, M_divergence=1.e3,
-            epsilon_stationnary=1.e-4, pc_threshold=0.1,
+            epsilon_stationnary=1.e-4,
             SA_params={"cutdeg": 8,
                        "SA_err": 1.e-6,
                        "cutdeg_glitch": 8,
diff --git a/examples/run_flake.py b/examples/run_flake.py
deleted file mode 100644
index adce1c70..00000000
--- a/examples/run_flake.py
+++ /dev/null
@@ -1,136 +0,0 @@
-# -*- coding: utf-8 -*-
-import numpy as np
-
-import fractalshades as fs
-import fractalshades.models as fsm
-import fractalshades.settings as settings
-import fractalshades.colors as fscolors
-
-def plot():
-    """
-    Example plot of "Dinkydau flake" location, classic test case for 
-    perturbation technique and glitch correction.
-    """
-    directory = "./flake"
-
-    # Dinkydau flake
-    # http://www.fractalforums.com/announcements-and-news/pertubation-theory-glitches-improvement/msg73027/#msg73027
-    # Ball method 1 found period: 7884
-    x = "-1.99996619445037030418434688506350579675531241540724851511761922944801584242342684381376129778868913812287046406560949864353810575744772166485672496092803920095332"
-    y = "0.00000000000000000000000000000000030013824367909383240724973039775924987346831190773335270174257280120474975614823581185647299288414075519224186504978181625478529"
-    dx = "1.8e-157"
-    precision = 200
-
-    # Set to True if you only want to rerun the post-processing part
-    settings.skip_calc = False
-    # Set to True to enable multi-processing
-    settings.enable_multiprocessing = True
-
-    nx = 3200
-    xy_ratio = 0.5
-    theta_deg = 0.
-    complex_type = np.complex128
-
-    mandelbrot = fsm.Perturbation_mandelbrot(directory)
-    mandelbrot.zoom(
-            precision=precision,
-            x=x,
-            y=y,
-            dx=dx,
-            nx=nx,
-            xy_ratio=xy_ratio,
-            theta_deg=theta_deg,
-            projection="cartesian",
-            antialiasing=False)
-
-    mandelbrot.calc_std_div(
-            complex_type=complex_type,
-            file_prefix="dev",
-            subset=None,
-            max_iter=50000,
-            M_divergence=1.e3,
-            epsilon_stationnary=1.e-3,
-            pc_threshold=0.1,
-            SA_params={"cutdeg": 8,
-                       "cutdeg_glitch": 8,
-                       "SA_err": 1.e-4,
-                       "use_Taylor_shift": True},
-            glitch_eps=1.e-6,
-            interior_detect=False,
-            glitch_max_attempt=20)
-
-    mandelbrot.run()
-
-    glitched = fs.Fractal_Data_array(mandelbrot, file_prefix="dev",
-                postproc_keys=('stop_reason', lambda x: x == 3), mode="r+raw")
-    potential_data_key = ("potential", {})
-
-
-    citrus2 = np.array([103, 189, 0]) / 255.
-    citrus_white = np.array([252, 251, 226]) / 255.
-
-    wheat1 = np.array([244, 235, 158]) / 255.
-    wheat2 = np.array([246, 207, 106]) / 255.
-    wheat3 = np.array([191, 156, 96]) / 255.
-
-    lavender1 = np.array([154, 121, 144]) / 255.
-    lavender2 = np.array([140, 94, 134]) / 255.
-    lavender3 = np.array([18, 16, 58]) / 255.
-    
-
-    def wave(x):
-        return 0.5 + (0.4 * (x - 0.5) - 0.6 * 0.5 * np.cos(x * np.pi * 3.))
-
-
-
-    colormap = fscolors.Fractal_colormap(
-        kinds="Lch",
-        colors1=np.vstack((citrus_white, wheat2, wheat1, wheat2, wheat1, wheat2, wheat3,
-                 wheat1, lavender2, wheat1, wheat2, wheat3, wheat1, lavender2, wheat1,
-                 lavender3, lavender2, lavender3, lavender1, lavender3, 
-                 lavender2)),
-        colors2=None,
-        n = 100,
-        funcs= lambda x: wave(x),
-        extent="mirror")
-        
-            
-
-
-#    colormap.extent = "mirror" #"repeat"
-
-    plotter = fs.Fractal_plotter(
-        fractal=mandelbrot,
-        base_data_key=potential_data_key, # potential_data_key, #glitch_sort_key,
-        base_data_prefix="dev",
-        base_data_function=lambda x:x,# np.sin(x*0.0001),
-        colormap=colormap,
-        probes_val=np.linspace(0., 1., 22) **0.2,# 200. + 200, #* 428  - 00.,#[0., 0.5, 1.], #phi * k * 2. + k * np.array([0., 1., 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) / 3.5,
-        probes_kind="qt",#"z", "qt"
-        mask=glitched)
-    
-    
-    #plotter.add_calculation_layer(postproc_key=potential_data_key)
-    
-    layer1_key = ("DEM_shade", {"kind": "potential",
-                                "theta_LS": 30.,
-                                "phi_LS": 50.,
-                                "shininess": 30.,
-                                "ratio_specular": 15000.})
-    plotter.add_grey_layer(postproc_key=layer1_key, intensity=0.75, 
-                         blur_ranges=[],#[[0.99, 0.999, 1.0]],
-                        disp_layer=False, #skewness=0.2,
-                         normalized=False, hardness=0.35,  
-            skewness=0.0, shade_type={"Lch": 1.0, "overlay": 1., "pegtop": 4.})
-    
-    layer2_key = ("field_lines", {})
-    plotter.add_grey_layer(postproc_key=layer2_key,
-                         hardness=1.0, intensity=0.68, skewness=0.4,
-                         blur_ranges=[[0.50, 0.60, 1.0]], 
-                         shade_type={"Lch": 0., "overlay": 2., "pegtop": 1.}) 
-
-
-    plotter.plot("dev", mask_color=(0., 0., 1.))
-
-if __name__ == "__main__":
-    plot()
diff --git a/examples/run_perturb2.py b/examples/run_perturb2.py
deleted file mode 100644
index 74b2c7e0..00000000
--- a/examples/run_perturb2.py
+++ /dev/null
@@ -1,131 +0,0 @@
-# -*- coding: utf-8 -*-
-import numpy as np
-
-import fractalshades as fs
-import fractalshades.models as fsm
-import fractalshades.settings as settings
-import fractalshades.colors as fscolors
-
-def plot():
-    """
-    Example plot of "Dinkydau flake" location, classic test case for 
-    perturbation techinque and glitch correction.
-    """
-    directory = "./perturb2"
-    # A simple showcas using perturbation technique
-    x, y = "-1.74920463345912691e+00", "-2.8684660237361114e-04"
-    dx = "5e-12"
-    precision = 16
-    nx = 800
-
-
-    # Set to True if you only want to rerun the post-processing part
-    settings.skip_calc = False
-    # Set to True to enable multi-processing
-    settings.enable_multiprocessing = True
-
-#    xy_ratio = 1.0
-#    theta_deg = 0.
-    # complex_type = np.complex128
-
-    mandelbrot = fsm.Perturbation_mandelbrot(directory)
-    mandelbrot.zoom(
-            precision=precision,
-            x=x,
-            y=y,
-            dx=dx,
-            nx=nx,
-            xy_ratio=1.0,
-            theta_deg=0., 
-            projection="cartesian",
-            antialiasing=False)
-
-    mandelbrot.calc_std_div(
-            complex_type=np.complex128,
-            file_prefix="dev",
-            subset=None,
-            max_iter=50000,
-            M_divergence=1.e3,
-            epsilon_stationnary=1.e-3,
-            pc_threshold=0.1,
-            SA_params={"cutdeg": 8,
-                       "cutdeg_glitch": 8,
-                       "SA_err": 1.e-4,
-                       "use_Taylor_shift": True},
-            glitch_eps=1.e-6,
-            interior_detect=True,
-            glitch_max_attempt=20)
-
-    mandelbrot.run()
-
-    cv = fs.Fractal_Data_array(mandelbrot, file_prefix="dev",
-                postproc_keys=('stop_reason', lambda x: x == 1), mode="r+raw")
-    potential_data_key = ("potential", {})
-
-
-
-
-
-    citrus2 = np.array([103, 189, 0]) / 255.
-    citrus_white = np.array([252, 251, 226]) / 255.
-
-    wheat1 = np.array([244, 235, 158]) / 255.
-    wheat2 = np.array([246, 207, 106]) / 255.
-    wheat3 = np.array([191, 156, 96]) / 255.
-
-    lavender1 = np.array([154, 121, 144]) / 255.
-    lavender2 = np.array([140, 94, 134]) / 255.
-    lavender3 = np.array([18, 16, 58]) / 255.
-    
-
-    def wave(x):
-        return 0.5 + (0.4 * (x - 0.5) - 0.6 * 0.5 * np.cos(x * np.pi * 3.))
-
-
-    colormap = fscolors.Fractal_colormap(
-        kinds="Lch",
-        colors1=np.vstack((citrus_white, wheat2, wheat1, wheat2, wheat1, wheat2, wheat3,
-                 wheat1, lavender2, wheat1, wheat2, wheat3, wheat1, lavender2, wheat1,
-                 lavender3, lavender2, lavender3, lavender1, lavender3, 
-                 lavender2)),
-        colors2=None,
-        n = 100,
-        funcs= lambda x: wave(x),
-        extent="mirror")
-
-
-    plotter = fs.Fractal_plotter(
-        fractal=mandelbrot,
-        base_data_key=potential_data_key, # potential_data_key, #glitch_sort_key,
-        base_data_prefix="dev",
-        base_data_function=lambda x:x,# np.sin(x*0.0001),
-        colormap=colormap,
-        probes_val=np.linspace(0., 1., 22) *0.5,# 200. + 200, #* 428  - 00.,#[0., 0.5, 1.], #phi * k * 2. + k * np.array([0., 1., 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) / 3.5,
-        probes_kind="qt",#"z", "qt"
-        mask=~cv)
-    
-    
-    #plotter.add_calculation_layer(postproc_key=potential_data_key)
-    
-    layer1_key = ("DEM_shade", {"kind": "potential",
-                                "theta_LS": 30.,
-                                "phi_LS": 50.,
-                                "shininess": 3.,
-                                "ratio_specular": 15000.})
-    plotter.add_grey_layer(postproc_key=layer1_key, intensity=0.75, 
-                         blur_ranges=[],#[[0.99, 0.999, 1.0]],
-                        disp_layer=False, #skewness=0.2,
-                         normalized=False, hardness=0.35,  
-            skewness=0.0, shade_type={"Lch": 1.0, "overlay": 1., "pegtop": 4.})
-    
-    layer2_key = ("field_lines", {})
-    plotter.add_grey_layer(postproc_key=layer2_key,
-                         hardness=1.0, intensity=0.68, skewness=0.4,
-                         blur_ranges=[[0.50, 0.60, 1.0]], 
-                         shade_type={"Lch": 0., "overlay": 2., "pegtop": 1.}) 
-
-
-    plotter.plot("dev", mask_color=(0., 0., 1.))
-
-if __name__ == "__main__":
-    plot()
diff --git a/install_from_local-static.sh b/install_from_local-static.sh
index 75bf56dc..001a487f 100755
--- a/install_from_local-static.sh
+++ b/install_from_local-static.sh
@@ -6,7 +6,7 @@
 # python3 -m pip install --prefix=/home/geoffroy/.local/lib/python3.8/site-packages --force-reinstall --no-deps --editable .
 
 # "Normal" installation
-python3 -m pip install --user --force-reinstall --no-deps ./dist/fractalshades_G_BILLOTEY-0.2.0-py3-none-any.whl
+python3 -m pip install --user --force-reinstall --no-deps ./dist/fractalshades_G_BILLOTEY-0.3b1-py3-none-any.whl
 
 # upload to test repository
 # python3 -m twine upload --repository testpypi dist/*
diff --git a/install_from_local.sh b/install_from_local.sh
deleted file mode 100755
index 55c3e9f8..00000000
--- a/install_from_local.sh
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-
-# Dev mode installation : --force reinstall -e
-# https://setuptools.readthedocs.io/en/latest/userguide/commands.html#develop
-# Note that --user option is not recognized for editable installs, hence the use of --prefix to explicitely specify the installation directory
-python3 -m pip install --prefix=/home/geoffroy/.local/lib/python3.8/site-packages --force-reinstall --no-deps --editable .
-
-# "Normal" installation
-# python3 -m pip install --user --force-reinstall --no-deps ./dist/fractalshades_G_BILLOTEY-0.2.0-py3-none-any.whl
-
-# upload to test repository
-# python3 -m twine upload --repository testpypi dist/*
-
diff --git a/src/fractalshades/__init__.py b/src/fractalshades/__init__.py
index e8b0968b..8f2c0f04 100644
--- a/src/fractalshades/__init__.py
+++ b/src/fractalshades/__init__.py
@@ -1,11 +1,7 @@
 # -*- coding: utf-8 -*-
-__version__ = "0.2.0"
-
-
-from .core import *
-from .perturbation import *
-
-#global imported_numba_xr
-#if not "imported_numba_xr" in globals():
-#    import fractalshades.numpy_utils.numba_xr as numba_xr
+__version__ = "0.3.b1"
 
+from .settings import get_figures, close
+from .utils import zoom_options, calc_options
+from .core import Fractal, Fractal_plotter, _Pillow_figure
+from .perturbation import PerturbationFractal
diff --git a/src/fractalshades/colors/color_mapping.py b/src/fractalshades/colors/color_mapping.py
index 97a8018e..a171e7d6 100644
--- a/src/fractalshades/colors/color_mapping.py
+++ b/src/fractalshades/colors/color_mapping.py
@@ -11,19 +11,6 @@
 
 import fractalshades.numpy_utils.expr_parser as fs_parser
 
-#import fractalshades.settings as fssettings
-#import fractalshades.utils as fsutils
-
-#def mkdir_p(path):
-#    """ Creates directory ; if exists does nothing """
-#    try:
-#        os.makedirs(path)
-#    except OSError as exc:
-#        if exc.errno == errno.EEXIST and os.path.isdir(path):
-#            pass
-#        else:
-#            raise exc
-
 
 class Color_tools():
     """ A bunch of staticmethods
@@ -44,7 +31,6 @@ class Color_tools():
     # tungsten-filament lighting with correlated color temperature of 2856 K. 
     A_ref_white = np.array([1.0985, 1.0000, 0.3558])
 
-
     @staticmethod
     def rgb_to_XYZ(rgb):
         arr = np.swapaxes(rgb, 0, 1)
@@ -368,16 +354,23 @@ def __init__(self, colors, kinds, grad_npts=32, grad_funcs="x",
         Fractal_colormap, concatenates (n-1) color gradients (passing
         through n colors).
 
+        Parameters
+        ==========
         colors : rgb float np.array of shape [n, 3]
+            The successives colors of the colormap
         kinds: arrays of string [n-1] , "Lab" or "Lch"
-        grad_npts : int np.array of size [n-1] : number of internal points
-            used, for each gradient
-        grad_funcs : [n-1] array of callables mapping [0, 1.] to [0., 1] passed
-            to each gradient. Callable passed as a string expression of x var
-            (ie, the func is the evaluation of "lambda x: " + expr).
+            The kind of gradient between colors n and n+1 : either linear in
+            Lab space of Lch space
+        grad_npts : int np.array of size [n-1]
+            number of internal points stores, gradient between colors n and n+1
+            a typical value is 32.
+        grad_funcs : [n-1] array of callables mapping [0, 1.] to [0., 1]
+            These are passed to each gradient. Callable passed as a string
+            expression of x var
+            (ie, the callable is the evaluation of "lambda x: " + expr).
             Default to identity.
-        extent : scalar, ["clip", "mirror", "repeat"] What to do with out of
-            range values.
+        extent : "clip" | "mirror" | "repeat"
+            What to do with out of range values.
         """      
         self.colors = colors = np.asarray(colors)
         self.n_grads = n_grads = colors.shape[0] - 1
@@ -401,22 +394,6 @@ def __init__(self, colors, kinds, grad_npts=32, grad_funcs="x",
             raise ValueError("Callable grad_funcs deprecated, use string")
 
 
-        # If color 2 not provided... we roll and cut the last item
-        # colors2 = np.roll(colors, -1, axis=0)
-#        colors1 = colors[:-1]
-#        colors2 = colors[1:]
-#        kinds = kinds[:-1]
-#        n = n[:-1]
-#        funcs = funcs[:-1]
-#        npts -= 1
-
-#        fc = sum(tuple(_Fractal_colormap(Color_tools.color_gradient(
-#                       kinds[ipt], colors[ipt, :], colors[ipt + 1, :], grad_npts[ipt],
-#                       grad_funcs[ipt])) for ipt in range(1, npts)),
-#                 start=_Fractal_colormap(Color_tools.color_gradient(
-#                       kinds[0], colors[0, :], colors[1, :], grad_npts[0],
-#                       grad_funcs[0])))
-
         # Should define 2 internal arrays
         self._n_interp_colors = sum(grad_npts) - n_grads + 1
         self._interp_colors = np.empty([self._n_interp_colors, 3],
@@ -444,19 +421,6 @@ def __init__(self, colors, kinds, grad_npts=32, grad_funcs="x",
 
 
 
-
-#    def params(self):
-#        return (self.kinds, self.colors1, self.colors2, self.funcs, 
-#                self.extent)
-    # Disables operations, as "constructor" vals will not follow
-#    def __neg__(self):
-#        raise NotImplementedError()
-#    def __add__(self, other):
-#        raise NotImplementedError()
-#    def __sub__(self, other):
-#        raise NotImplementedError()
-
-
     def __repr__(self):
         """ Return a string that can be evaluated to restaure the colormap
         """
@@ -588,11 +552,32 @@ def repeat(x, ext_min, ext_max):
         """
         return ext_min + ((x - ext_min) % (ext_max - ext_min))
 
+#: A colormap circling through black, blue, white, orange, black
+classic_colormap = Fractal_colormap(
+    colors=[[0.00784314, 0.01960784, 0.14509804],
+            [0.17647059, 0.10980392, 0.10588235],
+            [0.48627451, 0.24313725, 0.07058824],
+            [0.63921569, 0.39607843, 0.17647059],
+            [0.81176471, 0.58039216, 0.33333333],
+            [0.97647059, 0.85490196, 0.64313725],
+            [0.96470588, 0.98823529, 0.90196078],
+            [0.48627451, 0.7254902 , 0.90980392],
+            [0.27843137, 0.51764706, 0.74901961],
+            [0.12156863, 0.32941176, 0.57254902],
+            [0.07058824, 0.25490196, 0.49411765],
+            [0.00784314, 0.01960784, 0.14509804]],
+    kinds=['Lab', 'Lch', 'Lch', 'Lch', 'Lab', 'Lab', 'Lab', 'Lch', 'Lch', 'Lch', 'Lab'],
+    grad_npts=[32, 32, 32, 32,  32,  32,  32,  32,  32,  32,  32,  32],
+    grad_funcs=['x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x'],
+    extent='mirror'
+) 
+
 
 class Curve:
     def __init__(self, fn="x", brightness=None, hardness=None,
                  blur_ranges=None):
         """
+        TODO - work in progress
         A transfert curve from [0, 1] to [0, 1]
         ("Grey levels map")
         brightness float, O. leaves untouched, 0..+1 brighten
@@ -717,100 +702,3 @@ def blur(blur_ranges, data, blur_base):
             blur1, blur2, blur3 = blur_range
             data = data * Curve.bluring_coeff(blur1, blur2, blur3, blur_base)
         return data
-#
-#class _Fractal_colormap():
-#    # TODO : should be refactored so that it takes an array for which each line
-#    # is a list of argument to Color gradient array from Colortool class,
-#    # Each of them added to make the colormap
-#    # Discontinue support for matplotlib colormaps
-#    """
-#    Class responsible for mapping a real array to a colors array.
-#    Attributes :
-#        *_colors* Internal list of possible colors, [0 to self.n_colors]
-#        *_probes* list of indices in self._colors array, identifying the
-#                  transitions between differrent sections of the colormap. Each
-#                  of this probe is mapped to *_probe_value*, either given by the
-#                  user or computed at plot time.
-#    """
-#    def __init__(self, color_gradient, extent="mirror"):
-#        """
-#Creates a colormap from a color gradient array (as output by Color_tools
-#gradient functions, array of shape (n_colors, 3))
-#
-#*color_gradient*  a Color gradient array from Colortool class
-#*extent*  ["mirror", "repeat", "clip"] specifies what to do with out of range Values.
-#        """
-#        self._colors = color_gradient
-#        n_colors, _ = color_gradient.shape
-#        self._probes = np.array([0, n_colors-1], dtype=np.float32)
-##        self.quantiles_ref = None
-#        self.extent = extent
-#
-#    @property
-#    def n_colors(self):
-#        """ Total number of colors in the colormap. """
-#        return self._colors.shape[0]
-#
-#    @property
-#    def npts(self):
-#        return self.colors.shape[0]
-#
-#    @property
-#    def probes(self):
-#        """ Position of the "probes" ie transitions between the different parts
-#        of the colormap. Read-only. """
-#        return np.copy(self._probes)
-#
-#    def __neg__(self):
-#        """
-#        Returns a reversed colormap
-#        """
-#        other = _Fractal_colormap(self._colors[::-1, :])
-#        other._probes = self._probes[-1] - self._probes[::-1]
-#        return other
-#
-#    def __add__(self, other):
-#        """ Concatenates 2 Colormaps """
-#        fcm = _Fractal_colormap(np.vstack([self._colors, other._colors[1:]]))
-##        print("self._probes", self._probes)
-##        print("other._probes", other._probes)
-##        print("other._probes + (self._probes[-1])", other._probes + (self._probes[-1]))
-#        #self._probes[-1] += 0.5
-#        fcm._probes = np.concatenate([
-#            self._probes,
-#            (other._probes + (self._probes[-1]))[1:]
-#            ])
-#        return fcm
-#
-#    def __sub__(self, other):
-#        """ Sbstract a colormap ie adds its reversed version """
-#        return self.__add__(other.__neg__())
-
-
-
-## Defines the neutral element for cmap addition
-#EMPTY_CMAP = _Fractal_colormap(np.array([]).reshape([0, 3]))
-#EMPTY_CMAP._probes = np.array([-1])
-def test_print_cmap():
-    gold = np.array([255, 210, 66]) / 255.
-    black = np.array([0, 0, 0]) / 255.
-    purple = np.array([181, 40, 99]) / 255.
-    citrus2 = np.array([103, 189, 0]) / 255.
-    colors = np.vstack((citrus2[np.newaxis, :],
-                         purple[np.newaxis, :]))
-    colormap = Fractal_colormap(kinds="Lch", colors=colors,
-         grad_npts=200, grad_funcs="x", extent="mirror")
-    print(colormap)
-    
-    colors2 = np.vstack((citrus2[np.newaxis, :],
-                         black[np.newaxis, :],
-                         gold[np.newaxis, :],
-                         purple[np.newaxis, :]))
-    colormap2 = Fractal_colormap(kinds="Lch", colors=colors2,
-         grad_npts=[200, 20, 10], grad_funcs=["x", "x**6", "(1-x)"], extent="mirror")
-    print(colormap2)
-    print(repr(colormap2))
-
-
-if __name__ == "__main__":
-    test_print_cmap()
\ No newline at end of file
diff --git a/src/fractalshades/colors/layers.py b/src/fractalshades/colors/layers.py
index 27ccffdf..46fb24ad 100644
--- a/src/fractalshades/colors/layers.py
+++ b/src/fractalshades/colors/layers.py
@@ -46,11 +46,29 @@ class Virtual_layer:
     default_mask_color = None
 
     def __init__(self, postname, func, output=True):
-        """ Base class for all layer objects.
-        postname : str, key in linked plotter postproc dictionary
-        func : None | mapping applied to the post as prelim step
-        output : bool, if True : an image will be saved
-        """
+        """ 
+    Base class for all layer objects.
+    
+    Parameters
+    ----------
+    postname : str
+        Unique identifier key for this layer. When added to a 
+        `fractalshades.Fractal_plotter` plotter (see 
+        `fractalshades.Fractal_plotter.add_layer`) the layer can then be
+        accessed through this key :
+            
+            layer = plotter["postname"]
+        
+        The `postname` must match a `fractalshades.postproc.Postproc`
+        associated with this plotter.
+
+    func : None | callable
+        mapping applied during the post-processing initial step
+
+    output : bool
+        If True : an image will be saved when
+        `fractalshades.Fractal_plotter.add_layer` directive issued.
+    """
         self.postname = postname
         self.min = np.inf     # neutral value for np.nanmin
         self.max = -np.inf    # neutral value for np.nanmax
@@ -86,10 +104,17 @@ def dtype(self):
 
     def set_mask(self, layer, mask_color=None):
         """
+        Masks the base-layer
+
         Parameters
         ----------
-        layer : Layer_mask
-        mask_color : 3 or 4 - uplet or None
+        layer : `Bool_layer` | `Grey_layer`
+            The layer mask
+        mask_color : None | 3-uplet or 4_uplet of floats belonging to [0., 1.]
+            If None, a default value will be provided depending on the masked
+            array subclass
+            If a n-uplet is provided, it shall be with RGB or RGBA values for a
+            colored layer
         """
         if mask_color is None: # default class-implementation
             mask_color = self.default_mask_color
@@ -237,12 +262,27 @@ class Color_layer(Virtual_layer):
 
     def __init__(self, postname, func, colormap, probes_z=[0, 1],
                  probes_kind="relative", output=True):
-        """ A RGB or RGBA layer
-        func   a function or string
-        colormap a Colormap
-        probe_values a 2-array of floats
-        probes_kind "relative" / "absolute"
-        output boolean
+        """
+        A colored layer.
+        
+        Parameters
+        ----------
+        postname : str
+            passed to `Virtual_layer` constructor
+        func :
+            passed to `Virtual_layer` constructor
+        colormap : fractalshades.colors.Fractal_colormap
+            a colormap
+        probes_z : 2-floats list
+            The preprocessing affine rescaling. If [0., 1.] and probes_kind is 
+            `relative` (default) the minimum value of the field will be mapped
+            to 0. and the maximum to 1.
+        probes_kind : "relative" | "absolute"
+            The key for  probes_z values. If relative, they are expressed
+            relatively to min and max field values (0. <-> min, 1. <-> max) ;
+            if absolute they are used directly
+        output : bool
+            passed to `Virtual_layer` constructor
         """
         super().__init__(postname, func, output)
         # Will only become RGBA if 
@@ -271,32 +311,58 @@ def mode(self):
         else:
             raise ValueError(self.mask_kind)
 
-    def overlay(self, layer, blend_mode):
+    def overlay(self, layer, overlay_mode):
         """
-        public API
-        Layer : Color_layer / Greyscale_layer
+        Combine the output of current layer with an overlayed layer *after*
+        color-mapping.
+
+        Parameters
+        ----------
+        layer : `Color_layer` | `Grey_layer`
+            The layer to overlay. See `Overlay_mode` documentation for the
+            options and restrictions.
+        overlay_mode : `Overlay_mode`
+            Object describing the kind of overlay used.
         """
-        self._modifiers += [("overlay", layer, blend_mode)]
+        self._modifiers += [("overlay", layer, overlay_mode)]
 
     def shade(self, layer, lighting):
         """
-        public API
-        layer : Normal_map_layer
-        lighting : Lighting_model
+        Adds scene lighting to the layer through a normal map effect.
+
+        Parameters
+        ----------
+        layer : `Normal_map_layer`
+            The layer holding the normal vector field.
+        lighting : `Blinn_lighting`
+            Object describing the scene lighting.
         """
         self._modifiers += [("shade", layer, lighting)]
 
     def set_twin_field(self, layer, scale):
         """ 
-        public API
+        Combine two postprocessing in one.
 
-        Combine (add with scaling) the 2 arrays before color-mapping
-        
-        `layer` the combined layer (VirtualLayer)
-        `scale` the scaling coefficient ; if 1. both layer will have the same
-                min-max contribtion 
+        Combine the current layer with another layer *before* color-mapping.
+        (i.e. add directly the underlying float
+        `fractalshades.postproc.Postproc` arrays), tuning the intensity of
+        the effect with a scaling coefficient.
 
-        Note : The mask of the current layer will be applied to its twin 
+        Parameters
+        ----------
+        layer : `Virtual_Layer` instance
+            The combined layer
+        scale : float
+            The scaling coefficient ; if 1. both layer will have the same
+            min-max contribtion. If 0.1 the combined layer will only contribute
+            to 10 % (This is based on min and max values reached for each
+            layer)
+
+        Notes
+        -----
+
+        .. note::
+            Note : The mask of the current layer will be applied.
         """
         if layer is None:
             self._twin_field = None
@@ -409,8 +475,31 @@ class Grey_layer(Virtual_layer):
     default_mask_color = 0.
     mode = "L"
 
-    def __init__(self, postname, func, curve=None, probes_z=[0, 1],
+    def __init__(self, postname, func, curve=None, probes_z=[0., 1.],
                  probes_kind="relative", output=True):
+        """
+        A grey layer.
+        
+        Parameters
+        ----------
+        postname : str
+            passed to `Virtual_layer` constructor
+        func :
+            passed to `Virtual_layer` constructor
+        curve : None | callable
+            A mapping from [0, 1] to [0, 1] applied *after* rescaling (this is 
+            the equivalent of a colormap in greyscale)
+        probes_z : 2-floats list
+            The preprocessing affine rescaling. If [0., 1.] and probes_kind is 
+            `relative` (default) the minimum value of the field will be mapped
+            to 0. and the maximum to 1.
+        probes_kind : "relative" | "absolute"
+            The key for  probes_z values. If relative, they are expressed
+            relatively to min and max field values (0. <-> min, 1. <-> max) ;
+            if absolute they are used directly
+        output : bool
+            passed to `Virtual_layer` constructor
+        """
         super().__init__(postname, func, output)
         # Will only become RGBA if masked & mask_color has transparency
         self.curve = curve
@@ -419,10 +508,17 @@ def __init__(self, postname, func, curve=None, probes_z=[0, 1],
 
     def set_mask(self, layer, mask_color=None):
         """
+        Masks the base-layer
+
         Parameters
         ----------
-        layer : Layer_mask
-        mask_color : float
+        layer : Bool_layer | Grey_layer
+            The layer mask
+        mask_color : None | float belonging to [0., 1.]
+            If None, a default value will be provided depending on the masked
+            array subclass
+            If a value is provided, it shall be with :
+            a float belonging to [0. (black) ; 1. (white)] interval
         """
         if mask_color is None: # default class-implementation
             mask_color = self.default_mask_color
@@ -505,6 +601,13 @@ def apply_mask(self, crop, crop_size, mask_arr, mask_color):
             raise ValueError(self.mask_kind)
 
 class Disp_Layer(Grey_layer):
+    """
+    Grey layer with 32-bits precision.
+    
+    Same functionnality as a `Grey_layer` but the image export is with 32-bits
+    precision (useful when used as an height map for 3d processing like
+    Blender)
+    """
     mode = "I"
     def get_grey(self, arr):
         return np.int32(arr * 65535) # 2**16-1
@@ -519,7 +622,7 @@ def __init__(self, postname, max_slope=70, output=True):
         `max_slope`: maximal angle (in degree) of the normal map, will be used
                      for re-normalisation
         """
-        super().__init__(postname, None, output)
+        super().__init__(postname, None, colormap=None, output=output)
         self.max_slope = max_slope * np.pi / 180
     
     def crop(self, chunk_slice):
@@ -605,39 +708,54 @@ class Blinn_lighting:
     def __init__(self, k_Ambient, color_ambient):
         """
         This class holds the properties for the scene lightsources. Its
-        instances can be passed as a parameter to a Color_layer 'shade' method.
-        k_Ambient : ambient lighting pourcentage
-        color_ambient: a RGB 3-tuple, channels val in [0, 1]
+        instances can be passed as a parameter to `Color_layer.shade`.
+        
+        Parameters
+        ==========
+        k_Ambient : 
+            ambient lighting coefficient. A typical value is 0.2
+        color_ambient: 3-uplet float
+            a RGB 3-tuple, channels values in [0, 1] : the ambient lighting
+            color. Usually white (1., 1., 1.)
         """
         # ref : https://en.wikipedia.org/wiki/Blinn%E2%80%93Phong_reflection_model
         # http://www.codinglabs.net/article_physically_based_rendering_cook_torrance.aspx
         self.k_Ambient = k_Ambient
-        self.color_ambient = color_ambient
+        self.color_ambient = np.asarray(color_ambient)
         self.light_sources = []
 
     def add_light_source(self, k_diffuse, k_specular, shininess,
                          angles, coords=None, color=np.array([1., 1., 1.]),
                          material_specular_color=None):
         """
-        k_diffuse : diffuse lighting pourcentage, as RGB
-        k_specular : diffuse lighting pourcentage, as RGB
-        shininess : Phong exponent for specular model
-        angles : angles (theta, phi) of light source at infinity, in degree
-                      theta = 0 : aligned to x
-                      phi = 0 : normal to z (= in xy plane) 
-        coords : position of lightsource from image center,
-                     in pourcentage of current image
-                    (if None (default), source is at infinity ; if non-null,
-                    LS_angle has no effect)
+        Adds a lightsource to the scene
+        
+        Parameters
+        ==========
+        k_diffuse : 3-uplet float
+            diffuse lighting coefficients, as RGB
+        k_specular : 3-uplet float
+            diffuse lighting pourcentage, as RGB
+        shininess : float
+            Phong exponent for specular model
+        angles : 2-uplet float
+            angles (theta, phi) of light source at infinity, in degree
+            theta = 0 : aligned to x
+            phi = 0 : normal to z (= in xy plane) 
+        coords : 2-uplet float
+            position of lightsource from image center,
+            in pourcentage of current image
+            (if None (default), source is at infinity ; 
+            if non-null, LS_angle has no effect)
         """
         self.light_sources += [{
-            "k_diffuse": k_diffuse,
-            "k_specular": k_specular,
+            "k_diffuse": np.asarray(k_diffuse),
+            "k_specular": np.asarray(k_specular),
             "shininess": shininess,
             "_angles": angles, 
             "angles": tuple(a * np.pi / 180. for a in angles),
             "coords": coords,
-            "color": color,
+            "color": np.asarray(color),
             "material_specular_color": material_specular_color
         }]
 
@@ -690,19 +808,22 @@ class Overlay_mode:
     def __init__(self, mode, **mode_options):
         """
         This class holds the properties for overlaying. Its instances can be
-        passed as a parameter to a Color_layer 'overlay' method.
+        passed as a parameter to `Color_layer.overlay` method.
 
         Parameters
         ----------
         mode:   "alpha_composite" | "tint_or_shade"
                 the kind of overlay  applied
-                alpha_composite : alpha-compositing of 2 color layers (its is
-                the mask of the upper layer which will be considered for
-                compositing)
-                tint_or_shade : the color of the lower layer will be adjusted
-                (tinted of shaded) depending of the upper layer value (tinted
-                above 0.5, sahded below 0.5). The upper layer in this case
-                shall be a Grey_layer instance.
+                
+                - alpha_composite :
+                    alpha-compositing of 2 color layers (its is
+                    the mask of the upper layer which will be considered for
+                    compositing)
+                - tint_or_shade : 
+                    the color of the lower layer will be adjusted
+                    (tinted of shaded) depending of the upper layer value
+                    (tinted above 0.5, shaded below 0.5). The upper layer in
+                    this case shall be a `Grey_layer` instance.
 
         mode_options : a dict with the applicable option for each mode:
 
diff --git a/src/fractalshades/core.py b/src/fractalshades/core.py
index 5f039996..882a96fa 100644
--- a/src/fractalshades/core.py
+++ b/src/fractalshades/core.py
@@ -22,18 +22,43 @@
 
 
 
-"""
+class _Pillow_figure:
+    def __init__(self, img, pnginfo):
+        """
+        This class is a wrapper that can be used to redirect a Fractal_plotter
+        output, for instance when generating the documentation.
+        """
+        self.img = img
+        self.pnginfo = pnginfo
+
+    def save_png(self, im_path):
+        self.img.save(im_path,format="png", pnginfo=self.pnginfo)
 
-"""
 
-class Fractal_plotter():
+class Fractal_plotter:
     def  __init__(self, postproc_batch):
         """
-        A Fractal plotter
-        - point to a single Fractal
+        The base plotting class.
+        
+        A Fractal plotter :
+
+        - points to a single Fractal
         - can hold several postproc_batches each one point to a single
           (fracal, calculation) couple
-         - can hold several post-processing layers
+        - can hold several post-processing layers
+         
+        Parameters
+        ----------
+        postproc_batch
+            A single `fractalshades.postproc.Postproc_batch` or a list of 
+            these
+
+        Notes
+        -----
+
+        .. warning::
+            When passed multiples `Postproc_batches`, each one shall point to
+            the same `Fractal` object
         """
         # postproc_batchescan be single or an enumeration
         postproc_batches = postproc_batch
@@ -72,6 +97,7 @@ def add_postproc_batch(self, postproc_batch):
         Note : several postproc_batches are needed whenever different 
         calculations need to be combined in an output plot, as a postproc
         batch can only map to a unique calc_name.
+        :meta private:
         """
         if postproc_batch.fractal != self.fractal:
             raise ValueError("Attempt to add a postproc_batch from a different"
@@ -83,11 +109,23 @@ def add_postproc_batch(self, postproc_batch):
 
     def add_layer(self, layer):
         """
-        Adds a layer field
-        Layer postname shall have been already registered in one of the plotter
-        batches.        
-        When a layer is added, a link layer-> Fractal_plotter is 
-        created ; a layer can only point to a single plotter.
+        Adds a layer field to allow subsequent graphical operations or outputs.
+        
+        Parameters
+        ----------
+        layer : `fractalshades.colors.layers.Virtual_layer` or a derived class
+            The layer to add
+        
+        Notes
+        -----
+
+        .. warning::
+            Layer `postname` shall have been already registered in one of the
+            plotter batches.
+
+        .. warning::
+            When a layer is added, a link layer -> Fractal_plotter is 
+            created ; a layer can only point to a single `Fractal_plotter`.
         """
         postname = layer.postname
         if postname not in (list(self.postnames) + self.postnames_2d):
@@ -108,6 +146,12 @@ def __getitem__(self, layer_name):
                 layer_name, list(l.postname for l in self.layers)))
 
     def plot(self):
+        """
+        The base method to produce images.
+        
+        When called, it will got through all the instance-registered layers
+        and plot each layer for which the `output` attribute is set to `True`.
+        """
         self.store_postprocs()
         self.compute_scalings()
         self.write_postproc_report()
@@ -117,6 +161,7 @@ def plot(self):
 
     def store_postprocs(self):
         """ Computes and stores posprocessed data in a temporary mmap
+        :meta private:
         """
         self.open_temporary_mmap()
         inc_posproc_rank = 0
@@ -237,7 +282,10 @@ def save_tagged(self, img, img_path, tag_dict):
         pnginfo = PIL.PngImagePlugin.PngInfo()
         for k, v in tag_dict.items():
             pnginfo.add_text(k, str(v))
-        img.save(img_path, pnginfo=pnginfo)
+        if fssettings.output_mode == "Pillow":
+            fssettings.add_figure(_Pillow_figure(img, pnginfo))
+        else:
+            img.save(img_path, pnginfo=pnginfo)
 
     def push_layers_to_images(self):
         for i, layer in enumerate(self.layers):
@@ -259,7 +307,7 @@ def push_cropped(self, chunk_slice, layer, im):
 
 
 
-class Fractal():
+class Fractal:
     
     REPORT_ITEMS = [
         "chunk1d_begin",
@@ -277,22 +325,111 @@ class Fractal():
         "U",
         "stop_reason",
         "stop_iter"]
-    
-    """
-    Abstract base class
-    Derived classes are reponsible for computing, storing and reloading the raw
-    data composing a fractal image.
-    
-    ref:
-https://en.wikibooks.org/wiki/Pictures_of_Julia_and_Mandelbrot_Sets/The_Mandelbrot_set
-    """
 
     def __init__(self, directory):
         """
-        Parameters
-    *directory*   The working base directory
-    *complex_type*  numpy type or ("Xrange", numpy type)
-        """
+The base class for all escape-time fractals calculations.
+
+Derived class should implement the actual calculation methods used in the
+innner loop. This class provides the outer looping (calculation is run on 
+successive tiles), enables multiprocessing, and manage raw-result storing 
+and retrieving.
+
+Parameters
+----------
+directory : str
+    Path for the working base directory
+
+Attributes
+----------
+directory : str
+    the working directory
+iref : int
+    the reference point index, when using perturbation technique (if
+    not, shall be None)
+glitch_max_attempt : int
+    The maximal number of glitch correction loops
+subset
+    A boolean array-like of the size of the image, when False the
+    calculation is skipped for this point. It is usually the result
+    of a previous calculation that can be passed via a Fractal_array
+    wrapper.
+    If None (default) all points will be calculated.
+glitch_stop_index : int
+    A calculation can exit for several reasons which are tracked in 
+    `stop_reason` int-array (see below the description of the raw data
+    arrays). Values of `stop_reason` above `glitch_stop_index` indicate
+    glitched pixels, for which the calculation is deemed invalid.
+complex_type :
+    the datatype used for Z output arrays
+codes :
+    the string identifier codes for the saved arrays:
+    (`complex_codes`, `int_codes`, `termination_codes`)
+
+Notes
+-----
+
+These notes describe implementation details and should be useful mostly to
+advanced users when subclassing.
+
+.. note::
+
+    **Special methods**
+    
+    this class and its subclasses may define several methods decorated with
+    specific tags:
+    
+    @`fsutils.zoom_options`
+        decorates the methods used to define the zoom
+    @fsutils.calc_options
+        decorates the methods defining the calculation inner-loop
+    @fsutils.interactive 
+        decorates the methods that can be called
+        interactively from the GUI (right-click then context menu selection).
+        The coordinates of the click are passed to the called method.
+
+.. note::
+    
+    **Calculation parameters**
+
+    To lanch a calculation, call `~fractalshades.Fractal.run`. The parameters
+    from the last 
+    @ `fractalshades.zoom_options` call and last 
+    @ `fractalshades.calc_options` call will be used. 
+    They are stored as class attributes, above a list of such attributes and
+    their  meaning (non-exhaustive as derived class cmay define their own).
+    Note that they should normally not be directly accessed but defined in
+    derived class through zoom and calc methods.
+
+.. note::
+
+    **Saved data**
+    
+        The calculation results (raw output of the inner loop at exit) are
+        saved to disk and internally accessed during plotting phase through
+        memory-mapping. These are:
+
+        chunk_mask    
+            boolean - alias for `subset`
+            Saved to disk as ``calc_name``\_Z.arr in ``data`` folder
+        Z
+            Complex fields, several fields can be defined and accessed through
+            a field string identifier.
+            Saved to disk as ``calc_name``\_Z.arr in ``data`` folder
+        U
+            Integer fields, several fields can be defined and accessed through
+            a field string identifier.
+            Saved to disk as ``calc_name``\_U.arr in ``data`` folder
+        stop_reason
+            Byte codes: the reasons for loop exit (max iteration reached ?
+            overflow ? other ?) A string identifier
+            Saved to disk as ``calc_name``\_stop_reason.arr in ``data`` folder
+        stop_iter
+            Integer: iterations count at loop exit
+            Saved to disk as ``calc_name``\_stop_iter.arr in ``data`` folder
+
+        The string identifiers are stored in ``codes`` attributes.
+"""
         self.directory = directory
         self.iref = None # None when no reference point used / needed
         self.glitch_max_attempt = 0
@@ -320,23 +457,43 @@ def zoom(self, *,
              projection: str="cartesian",
              antialiasing: bool=False):
         """
+        Define and stores as class-attributes the zoom parameters for the next
+        calculation.
+        
         Parameters
-    *x* *y*  coordinates of the central point
-    *dx*     reference data range (on x axis). Definition depends of the
-        projection selected, for 'cartesian' (default) this is the total span.
-    *nx*     number of pixels of the image along x axis
-    *xy_ratio*  ratio of dx / dy and nx / ny
-    *theta_deg*    Pre-rotation of the calculation domain
-    *complex_type*  numpy type or ("Xrange", numpy type)
-    *projection*   "cartesian" "spherical" "exp_map"
+        ----------
+        x : float
+            x-coordinate of the central point
+        y : float 
+            y-coordinate of the central point
+        dx : float
+            span of the view rectangle along the x-axis
+        nx : int
+            number of pixels of the image along the x-axis
+        xy_ratio: float
+            ratio of dx / dy and nx / ny
+        theta_deg : float
+            Pre-rotation of the calculation domain, in degree
+        projection : "cartesian" | "spherical" | "exp_map"
+            Kind of projection used (default to cartesian)
+        antialiasing : bool
+            If True, some degree of randomization is applied
         """
         # We're all set, the job is done by `zoom_options` wrapper...
 
     def run(self):
         """
-        Lauch a full calculation
+        Lauch a full calculation.
+                
+        The parameters from the last 
+        @ `fractalshades.zoom_options`\-tagged method call and last
+        @ `fractalshades.calc_options`\-tagged method call will be used.
+
+        If calculation results are already there, the parameters will be
+        compared and if identical, the calculation will be skipped. This is
+        done for each tile and each glitch correction iteration, so i enables
+        calculation to restart from an unfinished status.
         """
-        
         if not(self.res_available()):
             # We write the param file and initialize the
             # memmaps for progress reports and calc arrays
@@ -373,12 +530,14 @@ def px(self):
     @property
     def multiprocess_dir(self):
         """ Directory used for multiprocess stdout stderr streams redirection
+        :meta private:
         """
         return os.path.join(self.directory, "multiproc_calc")
 
     @property
     def Xrange_complex_type(self):
         """ Return True if the data type is a xrange array
+        :meta private:
         """
         if type(self.complex_type) is tuple:
             type_modifier, _ = self.complex_type
@@ -402,6 +561,7 @@ def base_float_type(self):
     def params(self):
         """ Used to tag an output image or check if data is already computed
         and stored
+        :meta private:
         """
         software_params = {
                 "Software": "fractalshades " + fs.__version__,
@@ -422,7 +582,13 @@ def params(self):
 
 
     def clean_up(self, calc_name):
-        """ Deletes all data files associated with a given calc_name
+        """ Deletes all saved data files associated with a given `calc_name`.
+        
+        Parameters
+        ----------
+        calc_name : str
+            The string identifying the calculation run for which we want to
+            delete the files. 
         """
         for pattern in [
                 calc_name + "_*.arr",
@@ -442,7 +608,9 @@ def clean_up(self, calc_name):
     @property
     def pts_count(self):
         """ Return the total number of points for the current calculation 
-        taking into account the `subset` parameter """
+        taking into account the `subset` parameter
+        :meta private:
+        """
         if self.subset is not None:
 #            print("in pts_count", self.subset)
 #            print("in pts_count", np.sum(self.subset[None]),
diff --git a/src/fractalshades/gui/gui_cmap.py b/src/fractalshades/gui/gui_cmap.py
deleted file mode 100644
index 852229f9..00000000
--- a/src/fractalshades/gui/gui_cmap.py
+++ /dev/null
@@ -1,340 +0,0 @@
-# -*- coding: utf-8 -*-
-""" DEPRECATED """
-import inspect
-import typing
-import dataclasses
-
-
-from PyQt5 import QtCore
-from PyQt5.QtCore import Qt
-#from PyQt5.QtGui import QIcon
-from PyQt5.QtCore import pyqtSignal, pyqtSlot
-
-from PyQt5 import QtWidgets, QtGui
-from PyQt5.QtWidgets import (
-    QMainWindow,
-    QApplication,
-    QWidget,
-    QAction,
-    QHBoxLayout,
-    QVBoxLayout,
-    QCheckBox,
-    QLabel,
-    QMenuBar,
-    QToolBar,
-    QComboBox,
-    QGridLayout,
-    QSpacerItem,
-#    QSizePolicy,
-#    QFrame,
-    QSpinBox,
-    QGroupBox,
-    # QTableWidget,
-    QTableView,
-    QColorDialog,
-    QStyledItemDelegate,
-    QTableWidgetItem,
-    QStyleOptionViewItem
-)
-#from PyQt5.QtWidgets import (
-#    QComboBox,
-#    QSpinBox
-#)
-
-#from inspector import Func_inspector
-#import functools
-
-import fractalshades.colors as fscolors
-import numpy as np
-import functools
-import tempfile
-import os
-
-
-class ColorDelegate(QStyledItemDelegate):
-    def __init__(self, parent):
-        """ Custom cell delegate to display / edit a colors
-        parent : the QTableWidget
-        """
-        super().__init__(parent)
-
-    def createEditor(self, parent, option, index):
-        dialog = QColorDialog(None) #
-        dialog.setOption(QColorDialog.DontUseNativeDialog)
-        dialog.setCurrentColor(index.data(Qt.BackgroundRole))
-        dialog.setCustomColor(0, index.data(Qt.BackgroundRole))
-        # The returned editor widget should have Qt::StrongFocus
-        dialog.setFocusPolicy(Qt.StrongFocus)
-        dialog.setFocusProxy(parent)
-        return dialog
-
-    def setEditorData(self, editor, index):
-        color = index.data(Qt.BackgroundRole)
-        editor.setCurrentColor(color)
-
-    def setModelData(self, editor, model, index):
-        """ If modal QColorDialog result code is Accepted, save color"""
-        if editor.result():
-            color = editor.currentColor()
-            model.setData(index, color, Qt.BackgroundRole)
-
-    def paint(self, painter, option, index):
-        """ Fill with BackgroundRole color + red rectangle for selection."""
-        # After painting, you should ensure that the painter is returned to the
-        # state it was supplied in when this function was called
-        painter.save()
-        selected = bool(option.state & QtWidgets.QStyle.State_Selected)
-        print("color data", index.data(Qt.BackgroundRole))
-        color = index.data(Qt.BackgroundRole)
-        if color is None:
-            return
-        painter.fillRect(option.rect, index.data(Qt.BackgroundRole))
-        if selected:
-            rect = option.rect
-            rect.adjust(1, 1, -1, -1)
-            pen = QtGui.QPen(Qt.red)
-            pen.setWidth(2)
-            painter.setPen(pen)
-            painter.drawRect(rect)
-        painter.restore()
-
-
-class SpinDelegate(QStyledItemDelegate):
-    def __init__(self, parent, options):
-        """ Custom cell delegate to display / edit an int through combo box
-        parent : the QTableWidget
-        """
-        print("init SPIN", options)
-        super().__init__(parent)
-        self.min_val = options["min"]
-        self.max_val = options["max"]
-
-
-    def createEditor(self, parent, option, index):
-        print("c editor")
-        editor = QSpinBox(parent)
-        editor.setFrame(False)
-        editor.setMinimum(0)
-        editor.setMaximum(100)
-        return editor
-
-    def setEditorData(self, editor, index):
-        print("set editor data")
-        val = index.data(Qt.EditRole)
-        editor.setValue(val)
-
-    def setModelData(self, editor, model, index):
-        """ save int val to the model"""
-        print("set MODEL data")
-        val = editor.value()
-        print("set model data, val", val)
-        model.setData(index, val, Qt.EditRole)
-
-    def updateEditorGeometry(self, editor, option, index):
-        print("Update geom")
-        editor.setGeometry(option.rect)
-#    def paint(self, painter, option, index):
-#        """ Fill with BackgroundRole color + red rectangle for selection."""
-#        # After painting, you should ensure that the painter is returned to the
-#        # state it was supplied in when this function was called
-#        painter.save()
-#        selected = bool(option.state & QtWidgets.QStyle.State_Selected)
-#        painter.fillRect(option.rect, index.data(Qt.BackgroundRole))
-#        if selected:
-#            rect = option.rect
-#            rect.adjust(1, 1, -1, -1)
-#            pen = QtGui.QPen(Qt.red)
-#            pen.setWidth(2)
-#            painter.setPen(pen)
-#            painter.drawRect(rect)
-#        painter.restore()
-
-
-class Qcmap_editor(QWidget):
-    """
-    Widget of a cmap image 
-    """
-    def __init__(self, parent, cmap, minwidth=200, height=20):
-        super().__init__(parent)
-        self._cmap = cmap
-        
-        layout = QVBoxLayout()
-        layout.addWidget(self.add_param_box())
-        layout.addWidget(self.add_table_box())
-        # layout.addWidget(self.add_preview_box())
-        layout.addStretch(1)
-        self.setLayout(layout)
-        
-        self._wget_n.valueChanged.connect(functools.partial(
-                self.event_filter, "size"))
-        
-    def add_param_box(self):
-        
-        param_box = QGroupBox("Cmap parameters")
-        # Choice of number of lines
-        self._wget_n = QSpinBox(self)
-        self._wget_n.setRange(2, 256)
-        self._wget_n.setValue(len(self._cmap.colors))
-        # Choice of cmap "extent"
-        extent_choices = ["mirror", "repeat", "clip"]
-        self._wget_extent = QComboBox(self)
-        self._wget_extent.addItems(extent_choices)
-        # preview
-        self._preview = Qcmap_image(self, self._cmap)
-        param_layout = QHBoxLayout()
-        param_layout.addWidget(self._wget_n)#self._param_widget)
-        param_layout.addWidget(self._wget_extent)#self._param_widget)
-        param_layout.addWidget(self._preview, stretch=1)
-        param_box.setLayout(param_layout)
-        return param_box
-
-    def add_table_box(self):
-        table_box = QGroupBox("Cmap data")
-        table_layout = QHBoxLayout()
-
-        self._table = QTableView()# QTableWidget()
-        # COLUMNS : colors, kinds, n, funcs=None
-        # self._table.setColumnCount(4)
-        self._table.setStyleSheet('''
-                QTableView {
-                selection-background-color: white;
-                }
-                QTableView::item::selected {
-                  border: 2px solid red;
-                }
-            ''')
-        self.populate_table()
-
-        self._table.setItemDelegateForColumn(0, ColorDelegate(self._table))
-        self._table.setItemDelegateForColumn(2, SpinDelegate(self._table,
-                {"min": 0, "max":512}))
-#        self._table.setHorizontalHeaderLabels((
-#                "color",
-#                "kind",
-#                "grad_pts",
-#                "grad_func"))
-#        self._table.horizontalHeader().setSectionResizeMode(
-#                QtWidgets.QHeaderView.Stretch)
-        
-        table_layout = QHBoxLayout()
-        table_box.setLayout(table_layout)
-        table_layout.addWidget(self._table)#self._param_widget)
-        return table_box
-    
-    def populate_table(self):
-        n_rows = len(self._cmap.colors)
-        # self._table.setRowCount(n_rows)
-        model = QtGui.QStandardItemModel(n_rows, 2)
-        self._table.setModel(model)
-
-        for irow in range(n_rows):
-            # The "colors" col
-            val = self._cmap.colors[irow, :]
-            val = QtGui.QColor(*list(int(255 * f) for f in val))
-            index = model.index(irow, 0, QtCore.QModelIndex())
-            model.setData(index, val)
-            
-#            color_item = self._table.item(irow, 0)
-#            if color_item is None:
-#                color_item = QTableWidgetItem("")
-#            color_item.setData(Qt.BackgroundRole, QtGui.QColor(
-#                    *list(int(255 * f) for f in val)))
-#            self._table.setItem(irow, 0, color_item)
-
-        for irow in range(n_rows - 1): #
-            # The "grad_npts" col
-            val = self._cmap.grad_npts[irow]
-            print("###################val", val)
-            index = model.index(irow, 2, QtCore.QModelIndex())
-            model.setData(index, val)
-            print("################### SET DATA")
-            
-            
-            
-            
-
-    def event_filter(self, source, val):
-        print("event", source, val)
-
-
-class Qcmap_image(QWidget):
-    """
-    Widget of a cmap image with expanding width, fixed height
-    """
-    cmap_modified = pyqtSignal()
-
-    def __init__(self, parent, cmap, minwidth=200, height=20):
-        super().__init__(parent)
-        self._cmap = cmap
-        self.setMinimumWidth(minwidth)
-        self.setMinimumHeight(height)
-        self.setMaximumHeight(height)
-        self.setSizePolicy(QtWidgets.QSizePolicy.Minimum,
-                           QtWidgets.QSizePolicy.Expanding)
-
-    def paintEvent(self, evt):
-        size = self.size()
-        nx, ny = size.width(), size.height()
-        QtGui.QPainter(self).drawImage(0, 0, self._cmap.output_ImageQt(nx, ny))
-
-
-
-
-
-def getapp():
-    app = QtCore.QCoreApplication.instance()
-    if app is None:
-        app = QApplication([])
-    return app
-
-
-
-def test_cmap_widget():
-    
-    gold = np.array([255, 210, 66]) / 255.
-    black = np.array([0, 0, 0]) / 255.
-
-    colors = np.vstack((gold[np.newaxis, :],
-                         black[np.newaxis, :],
-                         gold[np.newaxis, :]))
-#    colors2 = np.vstack((black[np.newaxis, :],
-#                         gold[np.newaxis, :]))
-    kinds = ["Lch", "Lch"]
-    n = 100
-    funcs = ["x**6", "1.- (1. - x)**6"]
-
-    colormap = fscolors.Fractal_colormap(
-            colors, kinds, n, funcs, extent="clip")
-    print("probes", colormap._probes)
-    
-    
-    class Mywindow(QMainWindow):
-        def __init__(self):
-            super().__init__(parent=None)
-
-            self.setWindowTitle('Testing inspector')
-            tb = QToolBar(self)
-            self.addToolBar(tb)
-#            print_dict = QAction("print dict")
-            tb.addAction("print_dict")
-            
-            tb.actionTriggered[QAction].connect(self.on_tb_action)
-            #self.setWindowState(Qt.WindowMaximized)
-            # And don't forget to call setCentralWidget to your main layout widget.
-            icone =  Qcmap_editor(self, colormap) # fsgui.
-            self.setCentralWidget(icone)
-            self._icone = icone
-
-        def on_tb_action(self, qa):
-            print("ON ACTION qa", qa)
-            
-    app = getapp()
-    win = Mywindow()
-    win.show()
-    app.exec()
-    
-    
-    
-
-if __name__ == "__main__":
-    test_cmap_widget()
\ No newline at end of file
diff --git a/src/fractalshades/gui/guimodel.py b/src/fractalshades/gui/guimodel.py
index 68a34b7b..8a7ac0c9 100644
--- a/src/fractalshades/gui/guimodel.py
+++ b/src/fractalshades/gui/guimodel.py
@@ -974,7 +974,8 @@ def updateEditorGeometry(self, editor, option, index):
     def validate(self, index):
         val = self.modifier(index.data(Qt.DisplayRole))
         try:
-            return fs_parser.acceptable_expr(ast.parse(val, mode="eval"))
+            return fs_parser.acceptable_expr(ast.parse(val, mode="eval"),
+                                             safe_vars=["x"])
         except (SyntaxError):
             return False
 
@@ -998,8 +999,8 @@ def __init__(self, parent, cmap_presenter):
 
         layout = QVBoxLayout()
         layout.addWidget(self.add_param_box())
-        layout.addWidget(self.add_table_box())
-        layout.addStretch(1)
+        layout.addWidget(self.add_table_box(), stretch=1)
+#        layout.addStretch(1)
         self.setLayout(layout)
 
         self._wget_n.valueChanged.connect(functools.partial(
@@ -1045,7 +1046,7 @@ def populate_param_box(self):
 
     def add_table_box(self):
         table_box = QGroupBox("Cmap data")
-        table_layout = QHBoxLayout()
+        table_layout = QVBoxLayout()
 
         self._table = QTableWidget()
         # COLUMNS : colors, kinds, n, funcs=None
@@ -1068,15 +1069,20 @@ def add_table_box(self):
                 "kind",
                 "grad_pts",
                 "grad_func"))
-        self._table.horizontalHeader().setSectionResizeMode(
-                QtWidgets.QHeaderView.Stretch)
+        # Horizontal Headr: expand last column
+        h_header = self._table.horizontalHeader()
+        h_header.setSectionResizeMode(QtWidgets.QHeaderView.ResizeToContents)
+        h_header.setSectionResizeMode(3, QtWidgets.QHeaderView.Stretch)
+        # Vertical Header: ResizeToContents
+        self._table.verticalHeader().setSectionResizeMode(
+                QtWidgets.QHeaderView.ResizeToContents)
 
         # QAbstractItemView::SelectionMode QAbstractItemView::SingleSelection
         self._table.setSelectionMode(QtWidgets.QAbstractItemView.SingleSelection)
         
         table_layout = QHBoxLayout()
         table_box.setLayout(table_layout)
-        table_layout.addWidget(self._table)
+        table_layout.addWidget(self._table, stretch=1)
         return table_box
     
     def populate_table(self):
@@ -1314,6 +1320,7 @@ def __init__(self, parent, view_presenter): # im=None):#, xy_ratio=None):
 
         # Zoom rectangle disabled
         self._rect = None
+        self._rect_under = None
         self._drawing_rect = False
         self._dragging_rect = False
 
@@ -1392,7 +1399,7 @@ def reset_im(self):
 
 #        if self._qim is not None:
 #            self._group.removeFromGroup(self._qim)
-        for item in [self._qim, self._rect]:
+        for item in [self._qim, self._rect, self._rect_under]:
             if item is not None:
                 self._group.removeFromGroup(item)
 
@@ -1406,6 +1413,7 @@ def reset_im(self):
             self._qim = None
         
         self._rect = None
+        self._rect_under = None
         self._drawing_rect = False
 
     @staticmethod
@@ -1542,6 +1550,8 @@ def on_mouse_left_release(self, event):
             if (self._rect_pos0 == self._rect_pos1):
                 self._group.removeFromGroup(self._rect)
                 self._rect = None
+                self._group.removeFromGroup(self._rect_under)
+                self._rect_under = None
                 print("cancel drawing RECT")
                 self.cancel_drawing_rect()
             else:
@@ -1657,9 +1667,14 @@ def draw_rect(self, pos0, pos1):
         rectF = QtCore.QRectF(topleft, bottomRight)
         if self._rect is not None:
             self._rect.setRect(rectF)
+            self._rect_under.setRect(rectF)
         else:
+            self._rect_under = QGraphicsRectItem(rectF)
+            self._rect_under.setPen(QtGui.QPen(QtGui.QColor("red"), 0, Qt.SolidLine))
+            self._group.addToGroup(self._rect_under)
+
             self._rect = QGraphicsRectItem(rectF)
-            self._rect.setPen(QtGui.QPen(QtGui.QColor("red"), 0, Qt.DashLine))
+            self._rect.setPen(QtGui.QPen(QtGui.QColor("black"), 0, Qt.DotLine))
             self._group.addToGroup(self._rect)
 
     def selection_corners(self, pos0, pos1):
@@ -1668,8 +1683,8 @@ def selection_corners(self, pos0, pos1):
         diffy = abs(pos1.y() - pos0.y())
         # Enforce the correct ratio
         radius_sq = diffx ** 2 + diffy ** 2
-        diffx0 = math.sqrt(radius_sq / (1. + self.xy_ratio ** 2))
-        diffy0 = diffx0 * self.xy_ratio
+        diffx0 = math.sqrt(radius_sq / (1. + 1. / self.xy_ratio ** 2))
+        diffy0 = diffx0 / self.xy_ratio
         topleft = QtCore.QPointF(pos0.x() - diffx0, pos0.y() - diffy0)
         bottomRight = QtCore.QPointF(pos0.x() + diffx0, pos0.y() + diffy0)
         return topleft, bottomRight
diff --git a/src/fractalshades/models/mandelbrot.py b/src/fractalshades/models/mandelbrot.py
index 03637298..f8c784b0 100644
--- a/src/fractalshades/models/mandelbrot.py
+++ b/src/fractalshades/models/mandelbrot.py
@@ -7,12 +7,20 @@
 
 
 class Mandelbrot(fs.Fractal):
-    def __init__(self, *args, **kwargs):
+    def __init__(self, directory):
+        """
+        A standard power-2 Mandelbrot Fractal. 
+        
+        Parameters
+        ==========
+        directory : str
+            Path for the working base directory
+        """
+        super().__init__(directory)
         # default values used for postprocessing (potential)
         self.potential_kind = "infinity"
         self.potential_d = 2
         self.potential_a_d = 1.
-        super().__init__(*args, **kwargs)
 
     @fsutils.calc_options
     def base_calc(self, *,
@@ -23,13 +31,40 @@ def base_calc(self, *,
             epsilon_stationnary: float,
             datatype):
         """
-        Basic iterations for Mandelbrot standard set
+    Basic iterations for Mandelbrot standard set (power 2).
+
+    Parameters
+    ==========
+    calc_name : str
+         The string identifier for this calculation
+    subset : `fractalshades.postproc.Fractal_array`
+        A boolean array-like, where False no calculation is performed
+        If `None`, all points are calculated. Defaults to `None`.
+    max_iter : int
+        the maximum iteration number. If reached, the loop is exited with
+        exit code "max_iter".
+    M_divergence : float
+        The diverging radius. If reached, the loop is exited with exit code
+        "divergence"
+    epsilon_stationnary : float
+        A small float to early exit non-divergent cycles (based on
+        cumulated dzndz product). If reached, the loop is exited with exit
+        code "stationnary" (Those points should belong to Mandelbrot set
+        interior). A typical value is 1.e-3
+    datatype :
+        The dataype for operation on complex. Usually `np.complex128`
+        
+    Notes
+    =====
+    The following complex fields will be calculated: *zn* and its
+    derivatives (*dzndz*, *dzndc*, *d2zndc2*).
+    Exit codes are *max_iter*, *divergence*, *stationnary*.
         """
         complex_codes = ["zn", "dzndz", "dzndc", "d2zndc2"]
         int_codes = []
         stop_codes = ["max_iter", "divergence", "stationnary"]
         self.codes = (complex_codes, int_codes, stop_codes)
-        self.init_data_types(np.complex128)
+        self.init_data_types(datatype)
 
         self.potential_M = M_divergence
 
@@ -49,7 +84,7 @@ def iterate():
             @numba.njit
             def numba_impl(Z, U, c, stop_reason, n_iter):
                 if n_iter >= max_iter:
-                    stop_reason[0] = 0 # TODO should we return here
+                    stop_reason[0] = 0
 
                 Z[3] = 2 * (Z[3] * Z[0] + Z[2]**2)
                 Z[2] = 2 * Z[2] * Z[0] + 1.
@@ -77,24 +112,68 @@ def newton_calc(self, *,
             eps_newton_cv,
             datatype):
         """
-        Newton iterations for Mandelbrot standard set interior
+    Newton iterations for Mandelbrot standard set interior (power 2).
+
+    Parameters
+    ==========
+    calc_name : str
+         The string identifier for this calculation
+    subset : 
+        A boolean array-like, where False no calculation is performed
+        If `None`, all points are calculated. Defaults to `None`.
+    known_orders : None | int list 
+        If not `None`, only the integer listed of their multiples will be 
+        candidates for the cycle order, the other rwill be disregarded.
+        If `None` all order between 1 and `max_order` will be considered.
+    max_order : int
+        The maximum value for ccyle 
+    eps_newton_cv : float
+        A small float to qualify the convergence of the Newton iteration
+        Usually a fraction of a view pixel.
+    datatype :
+        The dataype for operation on complex. Usually `np.complex128`
         
-        References:
-          https://mathr.co.uk/blog/2013-04-01_interior_coordinates_in_the_mandelbrot_set.html
-          https://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
-          https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set_interior
-          https://fractalforums.org/fractal-mathematics-and-new-theories/28/interior-colorings/3352
+    Notes
+    =====
+    
+    .. note::
+
+        The following complex fields will be calculated:
+    
+        *zr*
+            A (any) point belonging to the attracting cycle
+        *attractivity*
+            The cycle attractivity (for a convergent cycle it is a complex
+            of norm < 1.)
+        *dzrdc*
+            Derivative of *zr*
+        *dattrdc*
+            Derivative of *attractivity*
+    
+        The following integer field will be calculated:
+    
+        *order*
+            The cycle order
+    
+        Exit codes are *max_order*, *order_confirmed*.
+
+    References
+    ==========
+    .. [1] <https://mathr.co.uk/blog/2013-04-01_interior_coordinates_in_the_mandelbrot_set.html>
+
+    .. [2] <https://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html>
         """
-        complex_codes = ["zr", "attractivity", "dzrdc", "dattrdc", "d2attrdc2"]
+        complex_codes = ["zr", "attractivity", "dzrdc", "dattrdc",
+                         "_partial", "_zn"]
         int_codes = ["order"]
         stop_codes = ["max_order", "order_confirmed"]
         self.codes = (complex_codes, int_codes, stop_codes)
-        self.init_data_types(np.complex128)
+        self.init_data_types(datatype)
 
         def initialize():
             def func(Z, U, c, chunk_slice):
-                print("#################### initialize")
-                Z[4, :] = 0.
+                Z[5, :] = 0.
+                Z[4, :] = 1.e6 # bigger than any reasonnable np.abs(zn)
                 Z[3, :] = 0.
                 Z[2, :] = 0.
                 Z[1, :] = 0.
@@ -110,8 +189,17 @@ def numba_impl(Z, U, c, stop_reason, n_iter):
                     stop_reason[0] = 0
                     return
 
-                # n is the candidate cycle order. Early exit if it is not a
-                # multiple of one of the known_orders
+                # If n is not a 'partial' for this point it cannot be the 
+                # cycle order : early exit
+                Z[5] = Z[5]**2 + c
+                m = np.abs(Z[5])
+                if m < Z[4].real:
+                    Z[4] = m # Cannot assign to the real part 
+                else:
+                    return
+
+                # Early exit if n it is not a multiple of one of the
+                # known_orders (provided by the user)
                 if known_orders is not None:
                     valid = False
                     for order in known_orders:
diff --git a/src/fractalshades/models/perturbation_mandelbrot.py b/src/fractalshades/models/perturbation_mandelbrot.py
index e6bf3167..553f1e3f 100644
--- a/src/fractalshades/models/perturbation_mandelbrot.py
+++ b/src/fractalshades/models/perturbation_mandelbrot.py
@@ -14,13 +14,20 @@
 
 class Perturbation_mandelbrot(fs.PerturbationFractal):
     
-    def __init__(self, *args, **kwargs):
+    def __init__(self, directory):
+        """
+        An arbitrary precision power-2 Mandelbrot Fractal. 
+
+        Parameters
+        ----------
+        directory : str
+            Path for the working base directory
+        """
+        super().__init__(directory)
         # Sets default values used for postprocessing (potential)
         self.potential_kind = "infinity"
         self.potential_d = 2
         self.potential_a_d = 1.
-        super().__init__(*args, **kwargs)
-
 
     @staticmethod
     def _ball_method1(c, px, maxiter, M_divergence):#, M_divergence):
@@ -283,26 +290,73 @@ def series_approx(self, SA_init, SA_loop, SA_params, iref, calc_name):
     @fsutils.calc_options
     def calc_std_div(self, *,
         calc_name: str,
-        complex_type,
+        datatype,
         subset,
         max_iter: int,
         M_divergence: float,
         epsilon_stationnary: float,
-        pc_threshold: float=0.1,
         SA_params=None,
         glitch_eps=None,
         interior_detect: bool=False,
         glitch_max_attempt: int=0):
         """
-        Computes the full data and derivatives
-        - "zn"
-        - "dzndz" (only if *interior_detect* is True)
-        - "dzndc"
+    Perturbation iterations (arbitrary precision) for Mandelbrot standard set
+    (power 2).
+
+    Parameters
+    ==========
+    calc_name : str
+         The string identifier for this calculation
+    datatype :
+        The dataype for operation on complex. Usually `np.complex128`
+    subset : 
+        A boolean array-like, where False no calculation is performed
+        If `None`, all points are calculated. Defaults to `None`.
+    max_iter : int
+        the maximum iteration number. If reached, the loop is exited with
+        exit code "max_iter".
+    M_divergence : float
+        The diverging radius. If reached, the loop is exited with exit code
+        "divergence"
+    epsilon_stationnary : float
+        EXPERIMENTAL for perturbation.
+        A small float to early exit non-divergent cycles (based on
+        cumulated dzndz product). If reached, the loop is exited with exit
+        code "stationnary" (Those points should belong to Mandelbrot set)
+        Used only if interior_detect is True
+    SA_params :
+        The dictionnary of parameters for Series-Approximation :
+
+        .. list-table:: 
+           :widths: 20 80
+           :header-rows: 1
+
+           * - keys
+             - values 
+           * - cutdeg
+             - int: polynomial degree used for first iteration  
+           * - cutdeg_glitch
+             - int: polynomial degree used for glitch correction 
+           * - SA_err
+             - float: maximal relative error before stopping SA 
+
+        if `None` SA is not activated.
+    glitch_eps : float
+        A small float to qualify the glitched pixel. Typical value 1.e-4
+    interior_detect : bool
+        EXPERIMENTAL for perturbation.
+        If True activates interior point detection
+    glitch_max_attempt: int
+        The maximal number of attempt for glitched pixel correction
+        
+    References
+    ==========
+    .. [1] <https://mathr.co.uk/blog/2021-05-14_deep_zoom_theory_and_practice.html>
+
+    .. [2] <http://www.fractalforums.com/announcements-and-news>
         
-        Note: if *interior_detect* is false we still allocate the arrays for
-        *dzndz* but we do not iterate.
         """
-        self.init_data_types(complex_type)
+        self.init_data_types(datatype)
         
         # used for potential post-processing
         self.potential_M = M_divergence
@@ -489,33 +543,31 @@ def numba_impl(Z, U, c, stop_reason, n_iter, SA_iter,
 
 
 
-
-
-
-
-
     @fsutils.calc_options
     def calc_fast(self, *,
         calc_name: str,
-        complex_type,
+        datatype,
         subset,
         max_iter: int,
         M_divergence: float,
         epsilon_stationnary: float,
-        pc_threshold: float=0.1,
         SA_params=None,
         glitch_eps=None,
         glitch_max_attempt: int=0):
+        
         """
-        Computes the full data and derivatives
-        - "zn"
-        - "dzndz" (only if *interior_detect* is True)
-        - "dzndc"
+        Faster options for perturbation iterations (arbitrary precision) of
+        Mandelbrot standard set (power 2).
         
-        Note: if *interior_detect* is false we still allocate the arrays for
-        *dzndz* but we do not iterate.
+        Refer to `calc_std_div` for the explanation of input parameters used.
+        Compared to `calc_std_div` this method should preferably be used for
+        exploration as it makes a few approximations and does not computes
+        derivatives
+        (Plotting options likes scene-lighting will not be available ;
+        field lines might show noticeable banding)
         """
-        self.init_data_types(complex_type)
+
+        self.init_data_types(datatype)
         
         # used for potential post-processing
         self.potential_M = M_divergence
@@ -689,15 +741,6 @@ def numba_impl(Z, U, c, stop_reason, n_iter, SA_iter,
 
 
 
-
-
-
-
-
-
-
-
-
 @numba.njit
 def SA_run(SA_loop, P0, n_iter, ref_path, kcX, SA_err_sq):
     SA_valid = True
diff --git a/src/fractalshades/perturbation.py b/src/fractalshades/perturbation.py
index 8bb25a35..1c2a4856 100644
--- a/src/fractalshades/perturbation.py
+++ b/src/fractalshades/perturbation.py
@@ -18,9 +18,20 @@ class PerturbationFractal(fs.Fractal):
 
     def __init__(self, directory):
         """
-        Sames args and kwargs as Fractal except :
-            - expect a "precision" integer : number of significant digits
-            - x and y input as strings
+The base class for escape-time fractals calculations implementing the 
+perturbation technique.
+
+Derived class should implement the actual calculation methods used in the
+innner loop, but also reference loop calculation, and optionnaly,
+Series approximation, glitch detection.
+This class provides the outer looping (calculation is run on 
+successive tiles), enables multiprocessing, and manage raw-result storing 
+and retrieving.
+
+Parameters
+----------
+directory : str
+    Path for the working base directory
         """
         super().__init__(directory)
 
@@ -35,6 +46,31 @@ def zoom(self, *,
              theta_deg: float,
              projection: str="cartesian",
              antialiasing: bool=False):
+        """
+        Define and stores as class-attributes the zoom parameters for the next
+        calculation.
+        
+        Parameters
+        ----------
+        precision : int
+            number of significant digits to use for reference point calculation
+        x : str or mpmath.mpf
+            x-coordinate of the central point
+        y : str or mpmath.mpf
+            y-coordinate of the central point
+        dx : str or mpmath.mpf
+            span of the view rectangle along the x-axis
+        nx : int
+            number of pixels of the image along the x-axis
+        xy_ratio: float
+            ratio of dx / dy and nx / ny
+        theta_deg : float
+            Pre-rotation of the calculation domain, in degree
+        projection : "cartesian"
+            Kind of projection used (only "cartesian" supported)
+        antialiasing : bool
+            If True, some degree of randomization is applied
+        """
         mpmath.mp.dps = precision
         # In case the user inputs were strings, we override with mpmath scalars
         self.x = mpmath.mpf(x)
@@ -235,22 +271,30 @@ def reload_SA(self, iref, calc_name):
             SA_params = pickle.load(tmpfile)
         return SA_params
 
-    def run(self):#, FP_loop, FP_params, SA_init, SA_loop, SA_params, max_iter,
-               # initialize, iterate, subset, codes,
-               # calc_name, pc_threshold, iref=0, glitch_stop_index=None,
-               # glitch_sort_key=None, glitch_max_attempt=0):
-        """
-        glitch_stop_index : All points with a stop_reason >= glitch_stop_index
-                are considered glitched (early ref exit glitch, or dynamic
-                glitch). Subset of pixel for which 
-                stop_reason == glitch_stop_index will be sorted according to 
-                minimum value of field `glitch_sort_key` the minimal  point
-                used as new reference.
-                if None, no glitch correction.
-                
-        glitch_sort_key : the complex (Z array) field where is stored our
-            'priority value' to select the next reference pixel.
+    def run(self):
         """
+        Launch a full perturbation calculation with glitch correction.
+
+        The parameters from the last 
+        @ `fractalshades.zoom_options`\-tagged method call and last
+        @ `fractalshades.calc_options`\-tagged method call will be used.
+
+        If calculation results are already there, the parameters will be
+        compared and if identical, the calculation will be skipped. This is
+        done for each tile and each glitch correction iteration, so i enables
+        calculation to restart from an unfinished status.
+        """
+        #  glitch_stop_index : All points with a stop_reason >= glitch_stop_index
+        #          are considered glitched (early ref exit glitch, or dynamic
+        #          glitch). Subset of pixel for which 
+        #          stop_reason == glitch_stop_index will be sorted according to 
+        #          minimum value of field `glitch_sort_key` the minimal  point
+        #          used as new reference.
+        #          if None, no glitch correction.
+        #                
+        #  glitch_sort_key : the complex (Z array) field where is stored our
+        #          'priority value' to select the next reference pixel.
+
         if not(self.res_available()):
             # We write the param file and initialize the
             # memmaps for progress reports and calc arrays
diff --git a/src/fractalshades/postproc.py b/src/fractalshades/postproc.py
index ca96c4d9..6b14c8b9 100644
--- a/src/fractalshades/postproc.py
+++ b/src/fractalshades/postproc.py
@@ -11,8 +11,22 @@
 
 class Postproc_batch:
     def __init__(self, fractal, calc_name):
-        """ Container for data common to different Postproc instances,
-        for a given chunk_slice"""
+        """
+        Container for several postprocessing objects (instances of `Postproc`)
+        which share the same Fractal & calculation (*calc_name*)
+        
+        It is usefull to have such container to avoid re-calculation of data
+        that can be shared between several post-proicessing objects.
+        
+        Parameters
+        ==========
+        fractal : fractalshades.Fractal instance
+            The shared fractal object
+        calc_name : str 
+            The string identifier for the shared calculation
+            
+        
+        """
         self.fractal = fractal
         self.calc_name = calc_name
         self.posts = dict()
@@ -24,18 +38,33 @@ def __init__(self, fractal, calc_name):
     def postnames(self):
         return self.posts.keys()
 
-    def add_postproc(self, name, postproc):
+    def add_postproc(self, postname, postproc):
         """
-        Add successive post and ensure they have consistent fractal & calc_name
+        Add successive postprocessing and ensure they have consistent
+        fractal & calc_name
+        
+        Parameters
+        ==========
+        postname : str
+            The string identifier that will be associated with this
+            post-processing object
+        postproc : `Postproc` instance
+            The post-processing object to be added
+
+        Notes
+        =====
+
+        .. warning::
+
+            A Postproc can only be added to one batch.
         """
-        if postproc.batch is not None: # and != self.fractal:
+        if postproc.batch is not None:
             raise ValueError("Postproc can only be linked to one batch")
-#        if postproc.calc_name is not None: #!= self.calc_name:
-#            raise ValueError("Postproc from a different calculation provided")
+
 
         postproc.link_batch(self) # Give access to batch members from postproc
         if postproc.field_count == 1:
-            self.posts[name] = postproc
+            self.posts[postname] = postproc
 
         elif postproc.field_count == 2:
             # We need to add 2 separate fields
@@ -45,9 +74,9 @@ def add_postproc(self, name, postproc):
             y_field.link_sibling(x_field)
             x_field.link_batch(self)
             y_field.link_batch(self)
-            self.posts[name + "_x"] = x_field
-            self.posts[name + "_y"] = y_field
-            self.postnames_2d += [name]
+            self.posts[postname + "_x"] = x_field
+            self.posts[postname + "_y"] = y_field
+            self.postnames_2d += [postname]
         else:
             raise ValueError("postproc.field_count bigger than 2: {}".format(
                     postproc.field_count))
@@ -74,15 +103,21 @@ def clear_chunk_data(self):
 
 
 class Postproc:
-    """ Abstract base class """
+    """
+    Abstract base class for all postprocessing objects.
+    """
     field_count = 1
-    def __init__(self): #, **post_dic):
+    def __init__(self):
         """ Filter None values : removing them from post_dic.
         """
         self.batch = None
 
     def link_batch(self, batch):
-        """ Link to Postproc group """
+        """
+        Link to Postproc group
+
+        :meta private:
+        """
         self.batch = batch
 
     @property
@@ -102,7 +137,10 @@ def context(self):
         return self.batch.context
 
     def ensure_context(self, key):
-        """ Check that the provided context contains the expected data
+        """
+        Check that the provided context contains the expected data
+
+        :meta private:
         """
         try:
             return self.context[key]
@@ -122,11 +160,17 @@ class Raw_pp(Postproc):
     def __init__(self, key, func=None):
         """
         A raw postproc provide direct access to the res data stored (memmap)
-        during calculation. (Similar to what does Fractal_array but within a
-        Postproc batch)
-        key : the raw data key
+        during calculation. (Similar to what does `Fractal_array`  but 
+        here within a `Postproc_batch`)
+
+        Parameters
+        ==========
+        key : 
+            The raw data key. Should be one of the *complex_codes*,
+            *int_codes*, *termination_codes* for this calculation.
         func: None | callable | a str of variable x (e.g. "np.sin(x)")
-              will be applied to the raw data if not None
+              will be applied as a pre-processing step to the raw data if not
+              `None`
         """
         super().__init__()
         self.key = key
@@ -166,16 +210,33 @@ def __getitem__(self, chunk_slice):
 class Continuous_iter_pp(Postproc):
     def __init__(self, kind=None, d=None, a_d=None, M=None, epsilon_cv=None):
         """
-        kind : "infinity" | "convergent" | "transcendent"
-        "d": degree of polynome d >= 2 if kind = "infinity"
-        "a_d": coeff of higher order monome if kind = "infinity"
-        "M": High number for sopping iteration if kind = "infinity"
-        "epsilon_cv": Small number for sopping iteration if kind = "convergent"
-        
-        Note: if parameters are not provided, will defaut to the attribute of 
-        associated fractal with name "potential_<attr>"
-        eg : self.fractal.potential_kind for "kind"
-        (see logic in potential_dic).
+        Return a continuous iteration number
+
+        Parameters
+        ==========
+        kind :  str "infinity" | "convergent" | "transcendent"
+            the number of orbit points used for the average
+        d : None | int
+            degree of polynome d >= 2 if kind = "infinity"
+        a_d : None | float
+            coeff of higher order monome if kind = "infinity"
+        M : None | float
+            High number corresponding to the criteria for stopping iteration
+            if kind = "infinity"
+        epsilon_cv : None | float
+            Small number corresponding to the criteria for stopping iteration
+            if kind = "convergent"
+
+        Notes
+        =====
+
+        .. note::
+
+            Usually it is not recommended to pass any parameter ; 
+            if a parameter <param>
+            is not provided, it will defaut to the attribute of 
+            the fractal with name "potential_<param>", which should be the 
+            recommended value.
         """
         super().__init__()
         post_dic = {
@@ -269,7 +330,30 @@ def __getitem__(self, chunk_slice):
 
 class Fieldlines_pp(Postproc):
     def __init__(self, n_iter=5, swirl=1., damping_ratio=0.1):
-        """ Field lines """
+        """
+        Return a continuous orbit-averaged angular value, allowing to reveal
+        fieldlines
+
+        Parameters
+        ==========
+        n_iter :  int
+            the number of orbit points used for the average
+        swirl : float between -1. and 1.
+            adds a random dephasing effect between each successive orbit point
+        damping_ratio : float
+            a geometric damping is used, damping ratio represent the ratio
+            between the scaling coefficient of the last orbit point and 
+            the first orbit point. 
+            For no damping (arithmetic mean) use 1.
+
+        Notes
+        =====
+
+        .. note::
+
+            The Continuous iteration number  shall have been calculated before
+            in the same `Postproc_batch`
+        """
         super().__init__()
         self.n_iter = n_iter
         self.swirl = swirl
@@ -329,7 +413,7 @@ def __getitem__(self, chunk_slice):
 
 class TIA_pp(Postproc):
     def __init__(self, n_iter=5):
-        """ Triangular average inequality """
+        """ Triangular average inequality - TODO """
         super().__init__()
         self.n_iter = n_iter
 
@@ -362,11 +446,37 @@ def __getitem__(self, chunk_slice):
 class DEM_normal_pp(Postproc):
     field_count = 2
     def __init__(self, kind="potential"):
+        
         """ 
-        This is a special case as it returns a complex array, hence 
-        field_count = 2. Internally 2 fields will be created (through helper
-        class XYCoord_wrapper_pp)
-        `kind`: "Milnor" | "potential"
+        Return the 2 components of the distance estimation method derivative
+        (*normal map*).
+
+        This postproc is normally used in combination with a
+        `fractalshades.colors.layers.Normal_map_layer`
+
+        Parameters
+        ==========
+        `kind`:  str
+            if "potential" (default) DEM is base on the potential.
+            For Mandelbrot power 2, "Milnor" option is also available which 
+            gives similar but esthetically slightly different results.
+
+        Notes
+        =====
+
+        .. note::
+
+            The Continuous iteration number  shall have been calculated before
+            in the same `Postproc_batch`
+
+            Alternatively, if kind="Milnor", the following raw complex fields
+            must be available from the calculation:
+
+                - "zn", "dzndc", "d2zndc2"
+
+        References
+        ==========
+        .. [1] <https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set>
         """
         super().__init__()
         self.kind = kind
@@ -454,10 +564,21 @@ class DEM_pp(Postproc):
     field_count = 1
     def __init__(self, px_snap=None):
         """ 
-        return the height map of the distance estimation method for the
-        potential associated with this postproc batch.
-        `px_snap`:  None | float   pixels at a distance from the fractal lower
-                    than `px_snap` (in pixel) will be snapped to 0.
+        Return a distance estimation of the pixel to the fractal set.
+        
+        Parameters
+        ==========
+        `px_snap`:  None | float
+            if not None, pixels at a distance from the fractal lower
+            than `px_snap` (expressed in image pixel) will be snapped to 0.
+
+        Notes
+        =====
+    
+        .. note::
+
+            The Continuous iteration number  shall have been calculated before
+            in the same `Postproc_batch`
         """
         super().__init__(px_snap)
 
@@ -491,10 +612,22 @@ class Attr_normal_pp(Postproc):
     field_count = 2
     def __init__(self):
         """ 
-        This is a special case as it returns a complex array, hence 
-        field_count = 2. Internally 2 fields will be created (through helper
-        class XYCoord_wrapper_pp)
-        Attr is the cycle attractivity
+        Return the 2 components of the cycle attractivity derivative (*normal
+        map*).
+
+        This postproc is normally used in combination with a
+        `fractalshades.colors.layers.Normal_map_layer`
+
+        Notes
+        =====
+    
+        .. note::
+
+            The following complex fields must be available from a previously
+            run calculation:
+
+                - "attractivity"
+                - "order" (optionnal)
         """
         super().__init__()
 
@@ -531,7 +664,24 @@ class Attr_pp(Postproc):
     field_count = 1
     def __init__(self, scale_by_order=False):
         """ 
-        Attr is the cycle attractivity
+        Return the cycle attractivity
+        
+        Parameters
+        ==========
+        scale_by_order : bool
+            If True return the attractivity divided by ccyle order (this
+            allows more realistic 3d-view exports)
+
+        Notes
+        =====
+    
+        .. note::
+
+            The following complex fields must be available from a previously
+            run calculation:
+
+                - "attractivity"
+                - "order" (optionnal)
         """
         super().__init__()
         self.scale_by_order = scale_by_order
@@ -543,27 +693,36 @@ def __getitem__(self, chunk_slice):
             termination_dic) = self.raw_data
         # Plotting the 'domed' height map for the cycle attractivity
         attr = Z[complex_dic["attractivity"], :]
-        r = U[int_dic["order"], :]
         abs2_attr = np.real(attr)**2 + np.imag(attr)**2
         abs2_attr = np.clip(abs2_attr, None, 1.)
         # divide by the cycle order
         val = np.sqrt(1. - abs2_attr) # / r
         if self.scale_by_order:
+            r = U[int_dic["order"], :]
             val *= (1. / r)
         return val, {}
 
 class Fractal_array:
     def __init__(self, fractal, calc_name, key, func=None):
         """
-        A Fractal_array provide direct access to the res data stored (memmap)
-        during calculation (Similar to what does Fractal_array but within a
-        Postproc batch)
-        fractal : The reference Fracal object
-        calc_name : the calculation name
-        key : the raw data key
+        Class which provide a direct access to the res data stored (memory
+        mapping)  during calculation.
+
+        One application, it can be passed as a *subset* parameter to most
+        calculation models, e.g. : `fractalshades.models.Mandelbrot.base_calc`
+
+        Parameters
+        ==========
+        fractal : `fractalshades.Fractal` derived class
+            The reference fracal object
+        calc_name : str 
+            The calculation name
+        key : 
+            The raw data key. Should be one of the *complex_codes*,
+            *int_codes*, *termination_codes* for this calculation.
         func: None | callable | a str of variable x (e.g. "np.sin(x)")
-              will be applied to the raw data if not None
-        A raw postproc 
+              will be applied as a pre-processing step to the raw data if not
+              `None`
         """
         self.fractal = fractal
         self.calc_name = calc_name
diff --git a/src/fractalshades/settings.py b/src/fractalshades/settings.py
index d8f997f4..19caad0d 100644
--- a/src/fractalshades/settings.py
+++ b/src/fractalshades/settings.py
@@ -2,6 +2,12 @@
 """
 General settings for the package
 """
+import sys
+
+# this is a pointer to the module object instance itself.
+this = sys.modules[__name__]
+
+
 enable_multiprocessing = True # True
 skip_calc = False
 
@@ -12,4 +18,27 @@
 chunk_size = 200
 
 # Number of bins for plots
-HIST_BINS = 100
\ No newline at end of file
+HIST_BINS = 100
+
+# Output mode : "png" | "Pillow"
+# if "Pillow" 
+output_mode = "png"
+
+# figures : module-level container used if output_mode == "Pillow"
+# -> the output images are stagged and not directly written to disk.
+figures = list()
+
+def add_figure(fig):
+    """ add a figure to the stagged images """
+    this.figures += [fig]
+
+def get_figures():
+    """ push the stagged images """
+    return this.figures
+
+def close(which):
+    """ delete the stagged images """
+    if which == "all":
+        this.figures = list()
+    else:
+        raise ValueError(which)
diff --git a/src/fractalshades/utils.py b/src/fractalshades/utils.py
index 1341e326..ce20f6cb 100644
--- a/src/fractalshades/utils.py
+++ b/src/fractalshades/utils.py
@@ -47,7 +47,10 @@ def wrapper(self, *args, **kwargs):
 
 
 class _store_kwargs_and_func_name:
-    """ Decorator for an instance method.
+    def __new__(cls, dic_name):
+        """ 
+    Decorator for an instance method.
+
     Parameters
     ----------
     dicname, str
@@ -70,7 +73,6 @@ class _store_kwargs_and_func_name:
         - this decorator will raise a ValueError if method is called with any
          positional arguments.
     """
-    def __new__(cls, dic_name):
         cls.dic_name = dic_name
         return object.__new__(cls)
 
@@ -94,9 +96,11 @@ def m_iter():
         return {name: method for name, method in m_iter()}
 
 
-# Decorator for the zooming method (normally only one per class)
+
+#: Decorator for the zooming method (normally only one per class)
 zoom_options = _store_kwargs("zoom_options")
 
-# Decorator for the calculation methods (can be several per class)
+#: Decorator for the calculation methods (can be several per class)
 calc_options = _store_kwargs_and_func_name("calc_options")
 
+


Test patch:
diff --git a/tests/test_layers.py b/tests/test_layers.py
index 6466d62d..66986084 100644
--- a/tests/test_layers.py
+++ b/tests/test_layers.py
@@ -369,7 +369,7 @@ def test_overlay1(self):
         # Overlay : alpha composite
         overlay_mode = Overlay_mode("alpha_composite")
         plotter[layer_name].shade(plotter["DEM_map"], light)
-        plotter[layer_name].overlay(plotter["attr"], blend_mode=overlay_mode)
+        plotter[layer_name].overlay(plotter["attr"], overlay_mode=overlay_mode)
 
         plotter.plot()
         self.layer = plotter[layer_name]
diff --git a/tests/test_perturbation.py b/tests/test_perturbation.py
index 038885b1..9afbd5f6 100644
--- a/tests/test_perturbation.py
+++ b/tests/test_perturbation.py
@@ -459,26 +459,24 @@ def calc(self, x, y, dx, precision, nx, complex_type, test_name, calc_name,
 
         if calc_fast:
             mandelbrot.calc_fast(
-                complex_type=complex_type,
+                datatype=complex_type,
                 calc_name=calc_name,
                 subset=None,
                 max_iter=max_iter,
                 M_divergence=1.e3,
                 epsilon_stationnary=1.e-3,
-                pc_threshold=0.1,
                 SA_params=SA_params,
                 glitch_eps=1.e-6,
                 glitch_max_attempt=glitch_max_attempt)
         else:
             mandelbrot.calc_std_div(
-                complex_type=complex_type,
+                datatype=complex_type,
                 calc_name=calc_name,
                 subset=None,
                 max_iter=max_iter,
                 M_divergence=1.e3,
                 epsilon_stationnary=1.e-3,
                 interior_detect=interior_detect,
-                pc_threshold=0.1,
                 SA_params=SA_params,
                 glitch_eps=1.e-6,
                 glitch_max_attempt=glitch_max_attempt)
@@ -502,7 +500,7 @@ def check_current_layer(self, err_max=0.01):
         self.assertTrue(err < err_max)
 
 if __name__ == "__main__":
-    full_test = False
+    full_test = True
     runner = unittest.TextTestRunner(verbosity=2)
     if full_test:
         runner.run(test_config.suite([Test_Perturbation_mandelbrot]))