diff --git a/docs/conf.py b/docs/conf.py
index c8e9bcbcd..ea18accd3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -18,6 +18,7 @@
 import shutil
 import re
 
+
 # -- Project information -----------------------------------------------------
 
 project = 'MRCPP'
@@ -49,12 +50,9 @@
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.
-#
-on_rtd = os.environ.get('READTHEDOCS') == 'True'
-if on_rtd:
-    html_theme = 'default'
-else:
-    html_theme = 'sphinx_rtd_theme'
+import sphinx_rtd_theme
+html_theme = "sphinx_rtd_theme"
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # The name of an image file (relative to this directory) to place at the top
 # of the sidebar.
diff --git a/docs/index.rst b/docs/index.rst
index f2f92ed5d..52832d79c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -28,7 +28,7 @@ The code can be found on `GitHub <https://github.com/MRChemSoft/mrcpp>`_.
    install.rst
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
    :caption: Application Program Interface
 
    mrcpp_api/introduction
@@ -42,6 +42,7 @@ The code can be found on `GitHub <https://github.com/MRChemSoft/mrcpp>`_.
 
 .. toctree::
    :maxdepth: 2
-   :caption: Programmers Manual
+   :caption: Programmer's Manual 
 
-   programmers_manual/clang_tidy
+   programmers_manual/overview
+   programmers_manual/clang_tidy.rst
diff --git a/docs/mrcpp_api/introduction.rst b/docs/mrcpp_api/introduction.rst
index 78b9d0a05..1dcaa1098 100644
--- a/docs/mrcpp_api/introduction.rst
+++ b/docs/mrcpp_api/introduction.rst
@@ -98,3 +98,5 @@ the capture list (square brackets), see e.g. `cppreference.com
 <http://en.cppreference.com/w/cpp/language/lambda>`_ for more
 details on lambda functions and how to use the capture list.
 
+
+
diff --git a/docs/mrcpp_api/mwoperators.rst b/docs/mrcpp_api/mwoperators.rst
index 3eb4cf972..ec8c691df 100644
--- a/docs/mrcpp_api/mwoperators.rst
+++ b/docs/mrcpp_api/mwoperators.rst
@@ -11,6 +11,8 @@ program by including:
 
     #include "MRCPP/MWOperators"
 
+.. doxygenclass:: mrcpp::MWOperator
+   :members:
 
 ConvolutionOperator
 -------------------
diff --git a/docs/programmers_manual/MWNode.rst b/docs/programmers_manual/MWNode.rst
new file mode 100644
index 000000000..de93dcecd
--- /dev/null
+++ b/docs/programmers_manual/MWNode.rst
@@ -0,0 +1,12 @@
+---------------------
+MWNode
+---------------------
+
+This is an introduction to the mwtree class. We write a small overarching summary of the class where we define the 
+algorithm/equation/structure reasoning for having this class or where it fits with the rest of the code.
+
+.. doxygenclass:: mrcpp::MWNode
+   :members:
+   :protected-members:
+   :private-members:
+
diff --git a/docs/programmers_manual/MWTree.rst b/docs/programmers_manual/MWTree.rst
new file mode 100644
index 000000000..5edc891c7
--- /dev/null
+++ b/docs/programmers_manual/MWTree.rst
@@ -0,0 +1,11 @@
+------------------
+MWTree 
+------------------
+
+This is an introduction to the mwtree class. We write a small overarching summary of the class where we define the 
+algorithm/equation/structure reasoning for having this class or where it fits with the rest of the code.
+
+.. doxygenclass:: mrcpp::MWTree
+   :members:
+   :protected-members:
+   :private-members:
\ No newline at end of file
diff --git a/docs/programmers_manual/clang_tidy.rst b/docs/programmers_manual/clang_tidy.rst
index ba68d2edc..bfe2e1000 100644
--- a/docs/programmers_manual/clang_tidy.rst
+++ b/docs/programmers_manual/clang_tidy.rst
@@ -3,10 +3,6 @@
    You can adapt this file completely to your liking, but it should at least
    contain the root `toctree` directive.
 
-==================================
-Programmers manual
-==================================
-
 ----------
 Clang-tidy
 ----------
diff --git a/docs/programmers_manual/overview.rst b/docs/programmers_manual/overview.rst
new file mode 100644
index 000000000..ff451b2aa
--- /dev/null
+++ b/docs/programmers_manual/overview.rst
@@ -0,0 +1,12 @@
+===============
+Manual overview
+===============
+The following sections show detailed documentation about the classes in MRCPP for programmers.
+
+TODO: maybe add some low level theory/figures/algorithms before showing classes, could help to understand some of the more abstract methods
+
+.. toctree::
+   :maxdepth: 1
+
+   MWTree
+   MWNode
diff --git a/src/operators/IdentityConvolution.cpp b/src/operators/IdentityConvolution.cpp
index eb0bd5278..8e9bea624 100644
--- a/src/operators/IdentityConvolution.cpp
+++ b/src/operators/IdentityConvolution.cpp
@@ -29,9 +29,10 @@
 
 namespace mrcpp {
 
-/** @returns New IdentityConvolution object
+/** @brief Constructor  of the IdentityConvolution object
+ *  @returns New IdentityConvolution object
  *  @param[in] mra: Which MRA the operator is defined
- *  @param[in] pr: Build precision, closeness to delta function
+ *  @param[in] prec: Build precision, closeness to delta function
  *  @details This will project a kernel of a single gaussian with
  *  exponent sqrt(10/build_prec).
  */
@@ -51,6 +52,19 @@ IdentityConvolution<D>::IdentityConvolution(const MultiResolutionAnalysis<D> &mr
     Printer::setPrintLevel(oldlevel);
 }
 
+/** @brief Constructor of the IdentityConvolution object in case of Periodic Boundary Conditions (PBC)
+ *  @returns New IdentityConvolution object
+ *  @param[in] mra: Which MRA the operator is defined
+ *  @param[in] prec: Build precision, closeness to delta function
+ *  @param[in] root: root scale of operator.
+ *  @param[in] reach: width at root scale (applies to periodic boundary conditions)
+ *  @details This will project a kernel of a single gaussian with
+ *           exponent sqrt(10/build_prec). This version of the constructor
+ *           is used for calculations within periodic boundary conditions (PBC). 
+ *           The \a root parameter is the coarsest negative scale at wich the operator
+ *           is applied. The \a reach parameter is the bandwidth of the operator at
+ *           the root scale. For details see \ref MWOperator
+ */
 template <int D>
 IdentityConvolution<D>::IdentityConvolution(const MultiResolutionAnalysis<D> &mra, double prec, int root, int reach)
         : ConvolutionOperator<D>(mra, root, reach) {
@@ -71,4 +85,4 @@ template class IdentityConvolution<1>;
 template class IdentityConvolution<2>;
 template class IdentityConvolution<3>;
 
-} // namespace mrcpp
\ No newline at end of file
+} // namespace mrcpp
diff --git a/src/operators/MWOperator.h b/src/operators/MWOperator.h
index 4907a13ca..4e3962fdc 100644
--- a/src/operators/MWOperator.h
+++ b/src/operators/MWOperator.h
@@ -32,6 +32,13 @@
 
 namespace mrcpp {
 
+/** @class MWOperator
+ *
+ * @brief Fixme
+ *
+ * @details Fixme
+ *
+ */
 template <int D>
 class MWOperator {
 public:
diff --git a/src/trees/MWNode.cpp b/src/trees/MWNode.cpp
index 5389da796..709966c58 100644
--- a/src/trees/MWNode.cpp
+++ b/src/trees/MWNode.cpp
@@ -40,8 +40,9 @@ using namespace Eigen;
 
 namespace mrcpp {
 
-/** MWNode default constructor.
- *  Should be used only by NodeAllocator to obtain
+/** @brief MWNode default constructor.
+ *
+ * @details Should be used only by NodeAllocator to obtain
  *  virtual table pointers for the derived classes. */
 template <int D>
 MWNode<D>::MWNode()
@@ -57,6 +58,13 @@ MWNode<D>::MWNode()
     MRCPP_INIT_OMP_LOCK();
 }
 
+/** @brief MWNode constructor.
+ *
+ * @param[in] tree: the MWTree the root node belongs to
+ * @param[in] idx: the NodeIndex defining scale and translation of the node
+ *
+ * @details Constructor for an empty node, given the corresponding MWTree and NodeIndex
+ */
 template <int D>
 MWNode<D>::MWNode(MWTree<D> *tree, const NodeIndex<D> &idx)
         : tree(tree)
@@ -70,6 +78,14 @@ MWNode<D>::MWNode(MWTree<D> *tree, const NodeIndex<D> &idx)
     MRCPP_INIT_OMP_LOCK();
 }
 
+/** @brief MWNode constructor.
+ *
+ * @param[in] tree: the MWTree the root node belongs to
+ * @param[in] rIdx: the integer specifying the corresponding root node
+ *
+ * @details Constructor for root nodes. It requires the corresponding
+ * MWTree and an integer to fetch the right NodeIndex
+ */
 template <int D>
 MWNode<D>::MWNode(MWTree<D> *tree, int rIdx)
         : tree(tree)
@@ -83,6 +99,14 @@ MWNode<D>::MWNode(MWTree<D> *tree, int rIdx)
     MRCPP_INIT_OMP_LOCK();
 }
 
+/** @brief MWNode constructor.
+ *
+ * @param[in] parent: parent node
+ * @param[in] cIdx: child index of the current node 
+ *
+ * @details Constructor for leaf nodes. It requires the corresponding
+ * parent and an integer to identify the correct child.
+ */
 template <int D>
 MWNode<D>::MWNode(MWNode<D> *parent, int cIdx)
         : tree(parent->tree)
@@ -96,8 +120,15 @@ MWNode<D>::MWNode(MWNode<D> *parent, int cIdx)
     MRCPP_INIT_OMP_LOCK();
 }
 
-/** MWNode copy constructor.
- *  Creates loose nodes and copy coefs */
+/** @brief MWNode copy constructor.
+ *
+ * @param[in] node: the original node
+ * @param[in] allocCoef: if true MW coefficients are allocated and copied from the original node
+ *
+ * @details Creates loose nodes and optionally copy coefs. The node
+ * does not "belong" to the tree: it cannot be accessed by traversing
+ * the tree.
+ */
 template <int D>
 MWNode<D>::MWNode(const MWNode<D> &node, bool allocCoef)
         : tree(node.tree)
@@ -127,18 +158,31 @@ MWNode<D>::MWNode(const MWNode<D> &node, bool allocCoef)
     MRCPP_INIT_OMP_LOCK();
 }
 
-/** MWNode destructor.
- * Recursive deallocation of a node and all its decendants */
+/** @brief MWNode destructor.
+ *
+ * @details Recursive deallocation of a node and all its decendants
+ */
 template <int D> MWNode<D>::~MWNode() {
     if (this->isLooseNode()) this->freeCoefs();
     MRCPP_DESTROY_OMP_LOCK();
 }
 
+/** @brief Dummy deallocation of MWNode coefficients.
+ *
+ * @details This is just to make sure this method never really gets
+ * called (derived classes must implement their own version). This was
+ * to avoid having pure virtual methods in the base class.
+ */
 template <int D> void MWNode<D>::dealloc() {
     NOT_REACHED_ABORT;
 }
 
-/** Allocate the coefs vector. Only used by loose nodes. */
+/** @brief Allocate the coefs vector.
+ *
+ * @details This is only used by loose nodes, because the loose nodes
+ * are not treated by the NodeAllocator class.
+ *
+ */
 template <int D> void MWNode<D>::allocCoefs(int n_blocks, int block_size) {
     if (this->n_coefs != 0) MSG_ABORT("n_coefs should be zero");
     if (this->isAllocated()) MSG_ABORT("Coefs already allocated");
@@ -151,7 +195,12 @@ template <int D> void MWNode<D>::allocCoefs(int n_blocks, int block_size) {
     this->setIsAllocated();
 }
 
-/** Deallocation of coefficients. Only used by loose nodes. */
+/** @brief Deallocate the coefs vector.
+ *
+ * @details This is only used by loose nodes, because the loose nodes
+ * are not treated by the NodeAllocator class.
+ *
+ */
 template <int D> void MWNode<D>::freeCoefs() {
     if (not this->isLooseNode()) MSG_ABORT("Only loose nodes here!");
 
@@ -164,6 +213,8 @@ template <int D> void MWNode<D>::freeCoefs() {
     this->clearIsAllocated();
 }
 
+/** @brief Printout of node coefficients
+ */
 template <int D> void MWNode<D>::printCoefs() const {
     if (not this->isAllocated()) MSG_ABORT("Node is not allocated");
     println(0, "\nMW coefs");
@@ -174,6 +225,8 @@ template <int D> void MWNode<D>::printCoefs() const {
     }
 }
 
+/** @brief wraps the MW coefficients into an eigen vector object
+ */
 template <int D> void MWNode<D>::getCoefs(Eigen::VectorXd &c) const {
     if (not this->isAllocated()) MSG_ABORT("Node is not allocated");
     if (not this->hasCoefs()) MSG_ABORT("Node has no coefs");
@@ -182,6 +235,9 @@ template <int D> void MWNode<D>::getCoefs(Eigen::VectorXd &c) const {
     c = VectorXd::Map(this->coefs, this->n_coefs);
 }
 
+/** @brief sets all MW coefficients and the norms to zero 
+ *
+ */
 template <int D> void MWNode<D>::zeroCoefs() {
     if (not this->isAllocated()) MSG_ABORT("Coefs not allocated " << *this);
 
@@ -190,27 +246,68 @@ template <int D> void MWNode<D>::zeroCoefs() {
     this->setHasCoefs();
 }
 
-/** Attach a set of coefs to this node. Only used locally (the tree is not aware of this). */
+/** @brief Attach a set of coefs to this node. Only used locally (the tree is not aware of this).
+ */
 template <int D> void MWNode<D>::attachCoefs(double *coefs) {
     this->coefs = coefs;
     this->setHasCoefs();
 }
 
+/** @brief assigns values to a block of coefficients
+ *
+ * @param[in] c: the input coefficients
+ * @param[in] block: the block index
+ * @param[in] block_size: size of the block
+ *
+ * @details a block is typically containing one kind of coefficients
+ * (given scaling/wavelet in each direction). Its size is then \f$
+ * (k+1)^D \f$ and the index is between 0 and \f$ 2^D-1 \f$.
+ */
 template <int D> void MWNode<D>::setCoefBlock(int block, int block_size, const double *c) {
     if (not this->isAllocated()) MSG_ABORT("Coefs not allocated");
     for (int i = 0; i < block_size; i++) { this->coefs[block * block_size + i] = c[i]; }
 }
 
+/** @brief adds values to a block of coefficients
+ *
+ * @param[in] c: the input coefficients
+ * @param[in] block: the block index
+ * @param[in] block_size: size of the block
+ *
+ * @details a block is typically containing one kind of coefficients
+ * (given scaling/wavelet in each direction). Its size is then \f$
+ * (k+1)^D \f$ and the index is between 0 and \f$ 2^D-1 \f$.
+ */
 template <int D> void MWNode<D>::addCoefBlock(int block, int block_size, const double *c) {
     if (not this->isAllocated()) MSG_ABORT("Coefs not allocated");
     for (int i = 0; i < block_size; i++) { this->coefs[block * block_size + i] += c[i]; }
 }
 
+/** @brief sets values of a block of coefficients to zero
+ *
+ * @param[in] block: the block index
+ * @param[in] block_size: size of the block
+ *
+ * @details a block is typically containing one kind of coefficients
+ * (given scaling/wavelet in each direction). Its size is then \f$
+ * (k+1)^D \f$ and the index is between 0 and \f$ 2^D-1 \f$.
+ */
 template <int D> void MWNode<D>::zeroCoefBlock(int block, int block_size) {
     if (not this->isAllocated()) MSG_ABORT("Coefs not allocated");
     for (int i = 0; i < block_size; i++) { this->coefs[block * block_size + i] = 0.0; }
 }
 
+/** @brief forward MW transform from this node to its children
+ *
+ * @param[in] overwrite: if true the coefficients of the children are
+ * overwritten. If false the values are summed to the already present
+ * ones.
+ *
+ * @details it performs forward MW transform inserting the result
+ * directly in the right place for each child node. The children must
+ * already be present and its memory allocated for this to work
+ * properly.
+ */
 template <int D> void MWNode<D>::giveChildrenCoefs(bool overwrite) {
     assert(this->isBranchNode());
     if (not this->isAllocated()) MSG_ABORT("Not allocated!");
@@ -236,6 +333,17 @@ template <int D> void MWNode<D>::giveChildrenCoefs(bool overwrite) {
     }
 }
 
+/** @brief forward MW transform to compute scaling coefficients of a single child
+ *
+ * @param[in] cIdx: child index
+ * @param[in] overwrite: if true the coefficients of the children are
+ * overwritten. If false the values are summed to the already present
+ * ones.
+ *
+ * @details it performs forward MW transform in place on a loose
+ * node. The scaling coefficients of the selected child are then
+ * copied/summed in the correct child node.
+ */
 template <int D> void MWNode<D>::giveChildCoefs(int cIdx, bool overwrite) {
 
     MWNode<D> node_i = *this;
@@ -258,6 +366,10 @@ template <int D> void MWNode<D>::giveChildCoefs(int cIdx, bool overwrite) {
 
 /** Takes a MWParent and generates coefficients, reverse operation from
  * giveChildrenCoefs */
+/** @brief backward MW transform to compute scaling/wavelet coefficients of a parent
+ *
+ * \warning This routine is only used in connection with Periodic Boundary Conditions
+ */
 template <int D> void MWNode<D>::giveParentCoefs(bool overwrite) {
     MWNode<D> node = *this;
     MWNode<D> &parent = getMWParent();
@@ -274,8 +386,12 @@ template <int D> void MWNode<D>::giveParentCoefs(bool overwrite) {
     parent.calcNorms();
 }
 
-/** Takes the scaling coefficients of the children and stores them consecutively
- * in the  given vector. */
+/** @brief Copy scaling coefficients from children to parent
+ * 
+ * @details Takes the scaling coefficients of the children and stores
+ * them consecutively in the corresponding block of the parent,
+ * following the usual bitwise notation.
+ */
 template <int D> void MWNode<D>::copyCoefsFromChildren() {
     int kp1_d = this->getKp1_d();
     int nChildren = this->getTDim();
@@ -286,6 +402,12 @@ template <int D> void MWNode<D>::copyCoefsFromChildren() {
     }
 }
 
+/** @brief Generates scaling cofficients of children
+ * 
+ * @details If the node is a leafNode, it takes the scaling&wavelet
+ * coefficients of the parent and it generates the scaling
+ * coefficients for the children
+ */
 template <int D> void MWNode<D>::threadSafeGenChildren() {
     MRCPP_SET_OMP_LOCK();
     if (isLeafNode()) {
@@ -295,15 +417,16 @@ template <int D> void MWNode<D>::threadSafeGenChildren() {
     MRCPP_UNSET_OMP_LOCK();
 }
 
-/** Coefficient-Value transform
+/** @brief Coefficient-Value transform
  *
- * This routine transforms the scaling coefficients of the node to the
+ * @details This routine transforms the scaling coefficients of the node to the
  * function values in the corresponding quadrature roots (of its children).
- * Input parameter = forward: coef->value.
- * Input parameter = backward: value->coef.
  *
- * NOTE: this routine assumes a 0/1 (scaling on children 0 and 1)
- *       representation, in oppose to s/d (scaling and wavelet). */
+ * @param[in] operation: forward (coef->value) or backward (value->coef).
+ *
+ * NOTE: this routine assumes a 0/1 (scaling on child 0 and 1)
+ *       representation, instead of s/d (scaling and wavelet).
+ */
 template <int D> void MWNode<D>::cvTransform(int operation) {
     int kp1 = this->getKp1();
     int kp1_dm1 = math_utils::ipow(kp1, D - 1);
@@ -392,10 +515,10 @@ void MWNode<D>::cvTransform(int operation) {
 }
 */
 
-/** Multiwavelet transform: fast version
+/** @brief Multiwavelet transform
   *
-  * Application of the filters on one node to pass from a 0/1 (scaling
-  * on children 0 and 1) representation to an s/d (scaling and
+  * @details Application of the filters on one node to pass from a 0/1 (scaling
+  * on child 0 and 1) representation to an s/d (scaling and
   * wavelet) representation. Bit manipulation is used in order to
   * determine the correct filters and whether to apply them or just
   * pass to the next couple of indexes. The starting coefficients are
@@ -408,9 +531,9 @@ void MWNode<D>::cvTransform(int operation) {
   * is formally faster than the other algorithm, the separation of the
   * three dimensions prevent the possibility to use the norm of the
   * operator in order to discard a priori negligible contributions.
-
-  * Luca Frediani, August 2006
-  * C++ version: Jonas Juselius, September 2009 */
+  *
+  *  * @param[in] operation: compression (s0,s1->s,d) or reconstruction (s,d->s0,s1).
+  */
 template <int D> void MWNode<D>::mwTransform(int operation) {
     int kp1 = this->getKp1();
     int kp1_dm1 = math_utils::ipow(kp1, D - 1);
@@ -450,19 +573,19 @@ template <int D> void MWNode<D>::mwTransform(int operation) {
     }
 }
 
-/** Set all norms to Undefined. */
+/** @brief Set all norms to Undefined. */
 template <int D> void MWNode<D>::clearNorms() {
     this->squareNorm = -1.0;
     for (int i = 0; i < this->getTDim(); i++) { this->componentNorms[i] = -1.0; }
 }
 
-/** Set all norms to zero. */
+/** @brief Set all norms to zero. */
 template <int D> void MWNode<D>::zeroNorms() {
     this->squareNorm = 0.0;
     for (int i = 0; i < this->getTDim(); i++) { this->componentNorms[i] = 0.0; }
 }
 
-/** Calculate and store square norm and component norms, if allocated. */
+/** @brief Calculate and store square norm and component norms, if allocated. */
 template <int D> void MWNode<D>::calcNorms() {
     this->squareNorm = 0.0;
     for (int i = 0; i < this->getTDim(); i++) {
@@ -472,7 +595,7 @@ template <int D> void MWNode<D>::calcNorms() {
     }
 }
 
-/** Calculate and return the squared scaling norm. */
+/** @brief Calculate and return the squared scaling norm. */
 template <int D> double MWNode<D>::getScalingNorm() const {
     double sNorm = this->getComponentNorm(0);
     if (sNorm >= 0.0) {
@@ -482,7 +605,7 @@ template <int D> double MWNode<D>::getScalingNorm() const {
     }
 }
 
-/** Calculate and return the squared wavelet norm. */
+/** @brief Calculate and return the squared wavelet norm. */
 template <int D> double MWNode<D>::getWaveletNorm() const {
     double wNorm = 0.0;
     for (int i = 1; i < this->getTDim(); i++) {
@@ -496,7 +619,7 @@ template <int D> double MWNode<D>::getWaveletNorm() const {
     return wNorm;
 }
 
-/** Calculate the norm of one component (NOT the squared norm!). */
+/** @brief Calculate the norm of one component (NOT the squared norm!). */
 template <int D> double MWNode<D>::calcComponentNorm(int i) const {
     if (this->isGenNode() and i != 0) return 0.0;
     assert(this->isAllocated());
@@ -515,9 +638,9 @@ template <int D> double MWNode<D>::calcComponentNorm(int i) const {
     return std::sqrt(sq_norm);
 }
 
-/** Update the coefficients of the node by a mw transform of the scaling
- * coefficients of the children. Option to overwrite or add up existing
- * coefficients. */
+/** @brief Update the coefficients of the node by a mw transform of the scaling
+ * coefficients of the children.
+ */
 template <int D> void MWNode<D>::reCompress() {
     if (this->isGenNode()) NOT_IMPLEMENTED_ABORT;
     if (this->isBranchNode()) {
@@ -529,8 +652,12 @@ template <int D> void MWNode<D>::reCompress() {
     }
 }
 
-/** Recurse down until an EndNode is found, and then crop children with
- * too high precision. */
+/** @brief Recurse down until an EndNode is found, and then crop children below the given precision threshold
+ *
+ * @param[in] prec: precision required
+ * @param[in] splitFac: factor used in the split check (larger factor means tighter threshold for finer nodes)
+ * @param[in] absPrec: flag to switch from relative (false) to absolute (true) precision.
+ */
 template <int D> bool MWNode<D>::crop(double prec, double splitFac, bool absPrec) {
     if (this->isEndNode()) {
         return true;
@@ -560,8 +687,11 @@ template <int D> void MWNode<D>::genParent() {
     NOT_REACHED_ABORT;
 }
 
-/** Recursive deallocation of children and all their descendants.
- * Leaves node as LeafNode and children[] as null pointer. */
+/** @brief Recursive deallocation of children and all their descendants.
+ *
+ * @details
+ * Leaves node as LeafNode and children[] as null pointer.
+ */
 template <int D> void MWNode<D>::deleteChildren() {
     if (this->isLeafNode()) return;
     for (int cIdx = 0; cIdx < getTDim(); cIdx++) {
@@ -576,7 +706,7 @@ template <int D> void MWNode<D>::deleteChildren() {
     this->setIsLeafNode();
 }
 
-/** Recursive deallocation of parent and all their forefathers. */
+/** @brief Recursive deallocation of parent and all their forefathers. */
 template <int D> void MWNode<D>::deleteParent() {
     if (this->parent == nullptr) return;
     MWNode<D> &parent = getMWParent();
@@ -586,6 +716,8 @@ template <int D> void MWNode<D>::deleteParent() {
     this->parent = nullptr;
 }
 
+
+/** @brief Deallocation of all generated nodes . */
 template <int D> void MWNode<D>::deleteGenerated() {
     if (this->isBranchNode()) {
         if (this->isEndNode()) {
@@ -596,6 +728,7 @@ template <int D> void MWNode<D>::deleteGenerated() {
     }
 }
 
+/** @brief returns the coordinates of the centre of the node */
 template <int D> Coord<D> MWNode<D>::getCenter() const {
     auto two_n = std::pow(2.0, -getScale());
     auto scaling_factor = getMWTree().getMRA().getWorldBox().getScalingFactors();
@@ -605,6 +738,7 @@ template <int D> Coord<D> MWNode<D>::getCenter() const {
     return r;
 }
 
+/** @brief returns the upper bounds of the D-interval defining the node  */
 template <int D> Coord<D> MWNode<D>::getUpperBounds() const {
     auto two_n = std::pow(2.0, -getScale());
     auto scaling_factor = getMWTree().getMRA().getWorldBox().getScalingFactors();
@@ -614,6 +748,7 @@ template <int D> Coord<D> MWNode<D>::getUpperBounds() const {
     return ub;
 }
 
+/** @brief returns the lower bounds of the D-interval defining the node  */
 template <int D> Coord<D> MWNode<D>::getLowerBounds() const {
     auto two_n = std::pow(2.0, -getScale());
     auto scaling_factor = getMWTree().getMRA().getWorldBox().getScalingFactors();
@@ -623,9 +758,11 @@ template <int D> Coord<D> MWNode<D>::getLowerBounds() const {
     return lb;
 }
 
-/** Routine to find the path along the tree.
+/** @brief Routine to find the path along the tree.
+ *
+ * @param[in] nIdx: the sought after node through its NodeIndex
  *
- * Given the translation indices at the final scale, computes the child m
+ * @details Given the translation indices at the final scale, computes the child m
  * to be followed at the current scale in oder to get to the requested
  * node at the final scale. The result is the index of the child needed.
  * The index is obtained by bit manipulation of of the translation indices. */
@@ -643,9 +780,11 @@ template <int D> int MWNode<D>::getChildIndex(const NodeIndex<D> &nIdx) const {
     return cIdx;
 }
 
-/** Routine to find the path along the tree.
+/** @brief Routine to find the path along the tree.
  *
- * Given a point in space, determines which child should be followed
+ * @param[in] r: the sought after node through the coordinates of a point in space
+ *
+ * @detailsGiven a point in space, determines which child should be followed
  * to get to the corresponding terminal node. */
 template <int D> int MWNode<D>::getChildIndex(const Coord<D> &r) const {
     assert(hasCoord(r));
@@ -660,6 +799,18 @@ template <int D> int MWNode<D>::getChildIndex(const Coord<D> &r) const {
     return cIdx;
 }
 
+/** @brief Returns the quadrature points in a given node
+ *
+ * @param[in,out] pts: quadrature points in a \f$ d \times (k+1) \f$ matrix form.
+ *
+ * @details The original quadrature points are fetched and then
+ * dilated and translated. For each cartesian direction \f$ \alpha =
+ * x,y,z... \f$ the set of quadrature points becomes \f$ x^\alpha_i =
+ * 2^{-n} (x_i + l^\alpha \f$. By taking all possible
+ * \f$(k+1)^d\f$ combinations, they will then define a d-dimensional
+ * grid of quadrature points.
+ *
+ */
 template <int D> void MWNode<D>::getPrimitiveQuadPts(MatrixXd &pts) const {
     int kp1 = this->getKp1();
     pts = MatrixXd::Zero(D, kp1);
@@ -672,6 +823,19 @@ template <int D> void MWNode<D>::getPrimitiveQuadPts(MatrixXd &pts) const {
     for (int d = 0; d < D; d++) pts.row(d) = sFac * (roots.array() + static_cast<double>(l[d]));
 }
 
+/** @brief Returns the quadrature points in a given node
+ *
+ * @param[in,out] pts: quadrature points in a \f$ d \times (k+1) \f$ matrix form.
+ *
+ * @details The original quadrature points are fetched and then
+ * dilated and translated to match the quadrature points in the
+ * children of the given node. For each cartesian direction \f$ \alpha = x,y,z... \f$
+ * the set of quadrature points becomes \f$ x^\alpha_i = 2^{-n-1} (x_i + 2 l^\alpha + t^\alpha) \f$, where \f$ t^\alpha =
+ * 0,1 \f$. By taking all possible \f$(k+1)^d\combinations \f$, they will
+ * then define a d-dimensional grid of quadrature points for the child
+ * nodes. 
+ *
+ */
 template <int D> void MWNode<D>::getPrimitiveChildPts(MatrixXd &pts) const {
     int kp1 = this->getKp1();
     pts = MatrixXd::Zero(D, 2 * kp1);
@@ -687,6 +851,16 @@ template <int D> void MWNode<D>::getPrimitiveChildPts(MatrixXd &pts) const {
     }
 }
 
+/** @brief Returns the quadrature points in a given node
+ *
+ * @param[in,out] pts: expanded quadrature points in a \f$ d \times
+ * (k+1)^d \f$ matrix form.
+ *
+ * @details The primitive quadrature points are used to obtain a
+ * tensor-product representation collecting all \f$ (k+1)^d \f$
+ * vectors of quadrature points.
+ *
+ */
 template <int D> void MWNode<D>::getExpandedQuadPts(Eigen::MatrixXd &pts) const {
     MatrixXd prim_pts;
     getPrimitiveQuadPts(prim_pts);
@@ -701,6 +875,16 @@ template <int D> void MWNode<D>::getExpandedQuadPts(Eigen::MatrixXd &pts) const
     if (D >= 4) NOT_IMPLEMENTED_ABORT;
 }
 
+/** @brief Returns the quadrature points in a given node
+ *
+ * @param[in,out] pts: expanded quadrature points in a \f$ d \times 
+ * 2^d(k+1)^d \f$ matrix form.
+ *
+ * @details The primitive quadrature points of the children are used to obtain a
+ * tensor-product representation collecting all \f$ 2^d (k+1)^d \f$
+ * vectors of quadrature points.
+ *
+ */
 template <int D> void MWNode<D>::getExpandedChildPts(MatrixXd &pts) const {
     MatrixXd prim_pts;
     getPrimitiveChildPts(prim_pts);
@@ -725,12 +909,16 @@ template <int D> void MWNode<D>::getExpandedChildPts(MatrixXd &pts) const {
     }
 }
 
-/** Const version of node retriever that NEVER generates.
+/** @brief Const version of node retriever that NEVER generates.
  *
+ * @param[in] idx: the requested NodeIndex
+ *
+ * @details
  * Recursive routine to find and return the node with a given NodeIndex.
  * This routine returns the appropriate Node, or a NULL pointer if
  * the node does not exist, or if it is a GenNode. Recursion starts at at this
- * node and ASSUMES the requested node is in fact decending from this node. */
+ * node and ASSUMES the requested node is in fact decending from this node.
+ */
 template <int D> const MWNode<D> *MWNode<D>::retrieveNodeNoGen(const NodeIndex<D> &idx) const {
     if (getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -745,12 +933,16 @@ template <int D> const MWNode<D> *MWNode<D>::retrieveNodeNoGen(const NodeIndex<D
     return this->children[cIdx]->retrieveNodeNoGen(idx);
 }
 
-/** Node retriever that NEVER generates.
+/** @brief Node retriever that NEVER generates.
+ *
+ * @param[in] idx: the requested NodeIndex
  *
+ * @details
  * Recursive routine to find and return the node with a given NodeIndex.
  * This routine returns the appropriate Node, or a NULL pointer if
  * the node does not exist, or if it is a GenNode. Recursion starts at at this
- * node and ASSUMES the requested node is in fact decending from this node. */
+ * node and ASSUMES the requested node is in fact decending from this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveNodeNoGen(const NodeIndex<D> &idx) {
     if (getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -765,6 +957,18 @@ template <int D> MWNode<D> *MWNode<D>::retrieveNodeNoGen(const NodeIndex<D> &idx
     return this->children[cIdx]->retrieveNodeNoGen(idx);
 }
 
+/** @brief Node retriever that returns requested Node or EndNode (const version).
+ *
+ * @param[in] r: the coordinates of a point in the node
+ * @param[in] depth: the depth which one needs to descend
+ *
+ * @details Recursive routine to find and return the node given the
+ * coordinates of a point in space.  This routine returns the
+ * appropriate Node, or the EndNode on the path to the requested node,
+ * and will never create or return GenNodes.  Recursion starts at at
+ * this node and ASSUMES the requested node is in fact decending from
+ * this node.
+ */
 template <int D> const MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const Coord<D> &r, int depth) const {
     if (getDepth() == depth or this->isEndNode()) { return this; }
     int cIdx = getChildIndex(r);
@@ -772,13 +976,18 @@ template <int D> const MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const Coord<D
     return this->children[cIdx]->retrieveNodeOrEndNode(r, depth);
 }
 
-/** Node retriever that return requested Node or EndNode.
+/** @brief Node retriever that returns requested Node or EndNode.
  *
- * Recursive routine to find and return the node with a given NodeIndex.
- * This routine returns the appropriate Node, or the EndNode on the
- * path to the requested node, and will never create or return GenNodes.
- * Recursion starts at at this node and ASSUMES the requested node is in fact
- * decending from this node. */
+ * @param[in] r: the coordinates of a point in the node
+ * @param[in] depth: the depth which one needs to descend
+ *
+ * @details Recursive routine to find and return the node given the
+ * coordinates of a point in space.  This routine returns the
+ * appropriate Node, or the EndNode on the path to the requested node,
+ * and will never create or return GenNodes.  Recursion starts at at
+ * this node and ASSUMES the requested node is in fact decending from
+ * this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const Coord<D> &r, int depth) {
     if (getDepth() == depth or this->isEndNode()) { return this; }
     int cIdx = getChildIndex(r);
@@ -786,6 +995,17 @@ template <int D> MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const Coord<D> &r,
     return this->children[cIdx]->retrieveNodeOrEndNode(r, depth);
 }
 
+/** @brief Node retriever that returns requested Node or EndNode (const version).
+ *
+ * @param[in] idx: the NodeIndex of the requested node
+ *
+ * @details Recursive routine to find and return the node given the
+ * coordinates of a point in space.  This routine returns the
+ * appropriate Node, or the EndNode on the path to the requested node,
+ * and will never create or return GenNodes.  Recursion starts at at
+ * this node and ASSUMES the requested node is in fact decending from
+ * this node.
+ */
 template <int D> const MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const NodeIndex<D> &idx) const {
     if (getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -800,6 +1020,18 @@ template <int D> const MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const NodeInd
     return this->children[cIdx]->retrieveNodeOrEndNode(idx);
 }
 
+/** @brief Node retriever that returns requested Node or EndNode.
+ *
+ * @param[in] idx: the NodeIndex of the requested node
+ *
+ * @details
+ * Recursive routine to find and return the node given the
+ * coordinates of a point in space.  This routine returns the
+ * appropriate Node, or the EndNode on the path to the requested node,
+ * and will never create or return GenNodes.  Recursion starts at at
+ * this node and ASSUMES the requested node is in fact decending from
+ * this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const NodeIndex<D> &idx) {
     if (getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -814,12 +1046,17 @@ template <int D> MWNode<D> *MWNode<D>::retrieveNodeOrEndNode(const NodeIndex<D>
     return this->children[cIdx]->retrieveNodeOrEndNode(idx);
 }
 
-/** Node retriever that ALWAYS returns the requested node.
+/** @brief Node retriever that ALWAYS returns the requested node.
+ *
+ * @param[in] r: the coordinates of a point in the node
+ * @param[in] depth: the depth which one needs to descend
  *
+ * @details
  * Recursive routine to find and return the node with a given NodeIndex.
  * This routine always returns the appropriate node, and will generate nodes
  * that does not exist. Recursion starts at this node and ASSUMES the
- * requested node is in fact decending from this node. */
+ * requested node is in fact decending from this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveNode(const Coord<D> &r, int depth) {
     if (depth < 0) MSG_ABORT("Invalid argument");
 
@@ -831,12 +1068,16 @@ template <int D> MWNode<D> *MWNode<D>::retrieveNode(const Coord<D> &r, int depth
     return this->children[cIdx]->retrieveNode(r, depth);
 }
 
-/** Node retriever that ALWAYS returns the requested node, possibly without coefs.
+/** @brief Node retriever that ALWAYS returns the requested node, possibly without coefs.
  *
+ * @param[in] idx: the NodeIndex of the requested node
+ *
+ * @details
  * Recursive routine to find and return the node with a given NodeIndex. This
  * routine always returns the appropriate node, and will generate nodes that
  * does not exist. Recursion starts at this node and ASSUMES the requested
- * node is in fact decending from this node. */
+ * node is in fact decending from this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveNode(const NodeIndex<D> &idx) {
     if (getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -854,10 +1095,14 @@ template <int D> MWNode<D> *MWNode<D>::retrieveNode(const NodeIndex<D> &idx) {
  *
  * WARNING: This routine is NOT thread safe! Must be used within omp critical.
  *
+ * @param[in] idx: the NodeIndex of the requested node
+ *
+ * @details
  * Recursive routine to find and return the node with a given NodeIndex. This
  * routine always returns the appropriate node, and will generate nodes that
  * does not exist. Recursion starts at this node and ASSUMES the requested
- * node is in fact related to this node. */
+ * node is in fact related to this node.
+ */
 template <int D> MWNode<D> *MWNode<D>::retrieveParent(const NodeIndex<D> &idx) {
     if (getScale() < idx.getScale()) MSG_ABORT("Scale error")
     if (getScale() == idx.getScale()) return this;
@@ -868,11 +1113,15 @@ template <int D> MWNode<D> *MWNode<D>::retrieveParent(const NodeIndex<D> &idx) {
     return this->parent->retrieveParent(idx);
 }
 
-/** Gives the norm (absolute value) of the node at the given NodeIndex.
+/** @brief Gives the norm (absolute value) of the node at the given NodeIndex.
  *
+ * @param[in] idx: the NodeIndex of the requested node
+ *
+ * @details
  * Recursive routine to find the node with a given NodeIndex. When an EndNode is
  * found, do not generate any new node, but rather give the value of the norm
- * assuming the function is uniformly distributed within the node. */
+ * assuming the function is uniformly distributed within the node.
+ */
 template <int D> double MWNode<D>::getNodeNorm(const NodeIndex<D> &idx) const {
     if (this->getScale() == idx.getScale()) { // we're done
         assert(getNodeIndex() == idx);
@@ -887,7 +1136,10 @@ template <int D> double MWNode<D>::getNodeNorm(const NodeIndex<D> &idx) const {
     return this->children[cIdx]->getNodeNorm(idx);
 }
 
-/** Test if a given coordinate is within the boundaries of the node. */
+/** @brief Test if a given coordinate is within the boundaries of the node.
+ *
+ * @param[in] r: point coordinates
+ */
 template <int D> bool MWNode<D>::hasCoord(const Coord<D> &r) const {
     double sFac = std::pow(2.0, -getScale());
     const NodeIndex<D> &l = getNodeIndex();
@@ -919,8 +1171,11 @@ template <int D> bool MWNode<D>::isCompatible(const MWNode<D> &node) {
     //    return true;
 }
 
-/** Test if the node is decending from a given NodeIndex, that is, if they have
- * overlapping support. */
+/** @brief Test if the node is decending from a given NodeIndex, that is, if they have
+ * overlapping support.
+ *
+ * @param[in] idx: the NodeIndex of the requested node
+ */
 template <int D> bool MWNode<D>::isAncestor(const NodeIndex<D> &idx) const {
     int relScale = idx.getScale() - getScale();
     if (relScale < 0) return false;
@@ -936,6 +1191,10 @@ template <int D> bool MWNode<D>::isDecendant(const NodeIndex<D> &idx) const {
     NOT_IMPLEMENTED_ABORT;
 }
 
+/** @brief printout ofm the node content.
+ *
+ * @param[in] o: the output stream
+ */
 template <int D> std::ostream &MWNode<D>::print(std::ostream &o) const {
     std::string flags = "       ";
     o << getNodeIndex();
@@ -963,7 +1222,8 @@ template <int D> std::ostream &MWNode<D>::print(std::ostream &o) const {
 
 /** @brief recursively set maxSquaredNorm and maxWSquareNorm of parent and descendants
  *
- * @details normalization is such that a constant function gives constant value,
+ * @details
+ * normalization is such that a constant function gives constant value,
  * i.e. *not* same normalization as a squareNorm
  */
 template <int D> void MWNode<D>::setMaxSquareNorm() {
@@ -980,7 +1240,8 @@ template <int D> void MWNode<D>::setMaxSquareNorm() {
         }
     }
 }
-/** @brief recursively reset maxSquaredNorm and maxWSquareNorm of parent and descendants to value -1 */
+/** @brief recursively reset maxSquaredNorm and maxWSquareNorm of parent and descendants to value -1
+ */
 template <int D> void MWNode<D>::resetMaxSquareNorm() {
     auto n = this->getScale();
     this->maxSquareNorm = -1.0;
diff --git a/src/trees/MWNode.h b/src/trees/MWNode.h
index e9901fde8..abc74147b 100644
--- a/src/trees/MWNode.h
+++ b/src/trees/MWNode.h
@@ -23,9 +23,6 @@
  * <https://mrcpp.readthedocs.io/>
  */
 
-/**
- *  Simple n-dimensional node
- */
 
 #pragma once
 
@@ -40,6 +37,20 @@
 
 namespace mrcpp {
 
+/** @class MWNode
+ *
+ * @brief Base class for Multiwavelet nodes
+ *
+ * @details A MWNode will contain the scaling and wavelet coefficients
+ * to represent functions or operators within a Multiwavelet
+ * framework. The nodes are in multidimensional. The dimensionality is
+ * set thoucgh the template parameter D=1,2,3. In addition to the
+ * coefficients the node contains metadata such as the scale, the
+ * translation index, the norm, pointers to parent node and child
+ * nodes, pointer to the corresponding MWTree etc... See memeber and
+ * data descriptions for details.
+ * 
+ */
 template <int D> class MWNode {
 public:
     MWNode(const MWNode<D> &node, bool allocCoef = true);
@@ -155,24 +166,24 @@ template <int D> class MWNode {
     friend class OperatorNode;
 
 protected:
-    MWTree<D> *tree{nullptr};
+    MWTree<D> *tree{nullptr};    ///< Tree the node belongs to
     MWNode<D> *parent{nullptr};  ///< Parent node
     MWNode<D> *children[1 << D]; ///< 2^D children
 
-    double squareNorm{-1.0};
-    double componentNorms[1 << D]; ///< 2^D components
-    double maxSquareNorm{-1.0};    // Largest squared norm among itself and descendants.
-    double maxWSquareNorm{-1.0};   // Largest wavelet squared norm among itself and descendants.
+    double squareNorm{-1.0};       ///< Squared norm of all 2^D (k+1)^D coefficients
+    double componentNorms[1 << D]; ///< Squared norms of the separeted 2^D components 
+    double maxSquareNorm{-1.0};    ///< Largest squared norm among itself and descendants.
+    double maxWSquareNorm{-1.0};   ///< Largest wavelet squared norm among itself and descendants.
                                    // NB: must be set before used.
-    double *coefs{nullptr};
+    double *coefs{nullptr};        ///< the 2^D (k+1)^D MW coefficients
     int n_coefs{0};
 
-    int serialIx{-1};       // index in serial Tree
-    int parentSerialIx{-1}; // index of parent in serial Tree, or -1 for roots
-    int childSerialIx{-1};  // index of first child in serial Tree, or -1 for leafnodes/endnodes
+    int serialIx{-1};       ///< index in serial Tree
+    int parentSerialIx{-1}; ///< index of parent in serial Tree, or -1 for roots
+    int childSerialIx{-1};  ///< index of first child in serial Tree, or -1 for leafnodes/endnodes
 
-    NodeIndex<D> nodeIndex;
-    HilbertPath<D> hilbertPath;
+    NodeIndex<D> nodeIndex;     ///< Scale and translation of the node
+    HilbertPath<D> hilbertPath; ///< To be documented
 
     MWNode();
     MWNode(MWTree<D> *tree, int rIdx);
diff --git a/src/trees/MWTree.cpp b/src/trees/MWTree.cpp
index 32f7f25bd..2a61836cd 100644
--- a/src/trees/MWTree.cpp
+++ b/src/trees/MWTree.cpp
@@ -23,12 +23,6 @@
  * <https://mrcpp.readthedocs.io/>
  */
 
-/**
- *  \date April 20, 2010
- *  CTCC, University of Troms
- *
- */
-
 #include "MWTree.h"
 
 #include "MWNode.h"
@@ -45,27 +39,39 @@ using namespace Eigen;
 
 namespace mrcpp {
 
-/** MWTree constructor with SerialTree storage for nodes.
- * Creates an empty tree object. Node construction and assignment of most of
- * the parameters are done in derived classes. */
+/** @brief MWTree constructor.
+ *
+ * @param[in] mra: the multiresolution analysis object
+ * @param[in] n: the name of the tree (only for printing purposes)
+ *
+ * @details Creates an empty tree object, containing only the set of
+ * root nodes. The information for the root node configuration to use
+ * is in the mra object which is passed to the constructor.
+ */
 template <int D>
 MWTree<D>::MWTree(const MultiResolutionAnalysis<D> &mra, const std::string &n)
         : MRA(mra)
-        , order(mra.getOrder())
-        , kp1_d(math_utils::ipow(mra.getOrder() + 1, D))
+        , order(mra.getOrder()) /// polynomial order
+        , kp1_d(math_utils::ipow(mra.getOrder() + 1, D)) ///nr of scaling coefficients \f$ (k+1)^D \f$
         , name(n)
         , squareNorm(-1.0)
         , rootBox(mra.getWorldBox()) {
     this->nodesAtDepth.push_back(0);
 }
 
-/** MWTree destructor. */
+/** @brief MWTree destructor. */
 template <int D> MWTree<D>::~MWTree() {
     this->endNodeTable.clear();
     if (this->nodesAtDepth.size() != 1) MSG_ERROR("Nodes at depth != 1 -> " << this->nodesAtDepth.size());
     if (this->nodesAtDepth[0] != 0) MSG_ERROR("Nodes at depth 0 != 0 -> " << this->nodesAtDepth[0]);
 }
 
+/** @brief Deletes all the nodes in the tree
+  *
+  * @details This method will recursively delete all the nodes,
+  * including the root nodes. Derived classes will call this method
+  * when the object is deleted.
+  */
 template <int D> void MWTree<D>::deleteRootNodes() {
     for (int i = 0; i < this->rootBox.size(); i++) {
         MWNode<D> &root = this->getRootMWNode(i);
@@ -77,10 +83,11 @@ template <int D> void MWTree<D>::deleteRootNodes() {
 
 /** @brief Remove all nodes in the tree
  *
- * @details Leaves the tree inn the same state as after construction, i.e.
- * undefined function containing only root nodes without coefficients.
- * The assigned memory (nodeChunks in NodeAllocator) is NOT released,
- * but is immediately available to the new function.
+ * @details Leaves the tree in the same state as after construction,
+ * i.e.  undefined tree structure containing only root nodes without
+ * coefficients.  The assigned memory, including branch and leaf
+ * nodes, (nodeChunks in NodeAllocator) is NOT released, but is
+ * immediately available to the new function.
  */
 template <int D> void MWTree<D>::clear() {
     for (int i = 0; i < this->rootBox.size(); i++) {
@@ -93,11 +100,11 @@ template <int D> void MWTree<D>::clear() {
     this->clearSquareNorm();
 }
 
-/** Calculate the squared norm of a function represented as a tree.
+/** @brief Calculate the squared norm \f$ ||f||^2_{\ldots} \f$ of a function represented as a tree.
  *
- * Norm is calculated using endNodes only, but if your endNodeTable is
- * incomplete (e.g. within growTree), the missing nodes must be given in the
- * input work vector. Involves an MPI reduction operation. */
+ * @details The norm is calculated using endNodes only. The specific
+ * type of norm which is computed will depend on the derived class
+ */
 template <int D> void MWTree<D>::calcSquareNorm() {
     double treeNorm = 0.0;
     for (int n = 0; n < this->getNEndNodes(); n++) {
@@ -108,6 +115,29 @@ template <int D> void MWTree<D>::calcSquareNorm() {
     this->squareNorm = treeNorm;
 }
 
+/** @brief Full Multiwavelet transform of the tree in either directions
+ *
+ * @param[in] type: TopDown (from roots to leaves) or BottomUp (from
+ * leaves to roots) which specifies the direction of the MW transform
+ * @param[in] overwrite: if true, the result will overwrite
+ * preexisting coefficients.
+ *
+ * @details It performs a Multiwavlet transform of the whole tree. The
+ * input parameters will specify the direction (upwards or downwards)
+ * and whether the result is added to the coefficients or it
+ * overwrites them. See the documentation for the #mwTransformUp 
+ * and #mwTransformDown for details.
+ * \f[ 
+ * \pmatrix{
+ * s_{nl}\\
+ * d_{nl}
+ * }
+ * \rightleftarrows \pmatrix{
+ * s_{n+1,2l}\\
+ * s_{n+1,2l+1}
+ * }
+ * \f]
+ */
 template <int D> void MWTree<D>::mwTransform(int type, bool overwrite) {
     switch (type) {
         case TopDown:
@@ -122,9 +152,15 @@ template <int D> void MWTree<D>::mwTransform(int type, bool overwrite) {
     }
 }
 
-/** Regenerate all s/d-coeffs by backtransformation, starting at the bottom and
- * thus purifying all coefficients. Option to overwrite or add up existing
- * coefficients of BranchNodes (can be used after operator application). */
+/** @brief Regenerates all s/d-coeffs by backtransformation
+ *
+ * @details It starts at the bottom of the tree (scaling coefficients
+ * of the leaf nodes) and it generates the scaling and wavelet
+ * coefficients if the parent node. It then proceeds recursively all the
+ * way up to the root nodes. This is generally used after a function
+ * projection to purify the coefficients obtained by quadrature at
+ * coarser scales which are therefore not precise enough.
+ */
 template <int D> void MWTree<D>::mwTransformUp() {
     std::vector<MWNodeVector<D>> nodeTable;
     tree_utils::make_node_table(*this, nodeTable);
@@ -142,9 +178,17 @@ template <int D> void MWTree<D>::mwTransformUp() {
     }
 }
 
-/** Regenerate all scaling coeffs by MW transformation of existing s/w-coeffs
- * on coarser scales, starting at the rootNodes. Option to overwrite or add up
- * existing scaling coefficients (can be used after operator application). */
+/** @brief Regenerates all scaling coeffs by MW transformation of existing s/w-coeffs
+ * on coarser scales
+ *
+ * @param[in] overwrite: if true the preexisting coefficients are overwritten
+ *
+ * @details The transformation starts at the rootNodes and proceeds
+ * recursively all the way to the leaf nodes. The existing scaling
+ * coefficeints will either be overwritten or added to. The latter
+ * operation is generally used after the operator application.
+ *
+ */
 template <int D> void MWTree<D>::mwTransformDown(bool overwrite) {
     std::vector<MWNodeVector<D>> nodeTable;
     tree_utils::make_node_table(*this, nodeTable);
@@ -169,9 +213,12 @@ template <int D> void MWTree<D>::mwTransformDown(bool overwrite) {
     }
 }
 
-/** @brief Set the MW coefficients to zero, fixed grid
- * @details Keeps the node structure of the tree, even though the zero function
- * is representable at depth zero. Use cropTree to remove unnecessary nodes.*/
+/** @brief Set the MW coefficients to zero, keeping the same tree structure
+ *   
+ * @details Keeps the node structure of the tree, even though the zero
+ * function is representable at depth zero. One should then use \ref cropTree to remove
+ * unnecessary nodes.
+ */
 template <int D> void MWTree<D>::setZero() {
     TreeIterator<D> it(*this);
     while (it.next()) {
@@ -181,7 +228,10 @@ template <int D> void MWTree<D>::setZero() {
     this->squareNorm = 0.0;
 }
 
-/** Increment node counters for non-GenNodes. This routine is not thread
+/** @brief Increments node counter by one for non-GenNodes.
+ *
+ * @details TO BE DOCUMENTED
+ * \warning: This routine is not thread
  * safe, and must NEVER be called outside a critical region in parallel.
  * It's way. way too expensive to lock the tree, so don't even think
  * about it. */
@@ -202,10 +252,14 @@ template <int D> void MWTree<D>::incrementNodeCount(int scale) {
     }
 }
 
-/** Decrement node counters for non-GenNodes. This routine is not thread
+/** @brief Decrements node counter by one for non-GenNodes.
+ *
+ * @details TO BE DOCUMENTED
+ * \warning: This routine is not thread
  * safe, and must NEVER be called outside a critical region in parallel.
  * It's way. way too expensive to lock the tree, so don't even think
- * about it. */
+ * about it.
+ */
 template <int D> void MWTree<D>::decrementNodeCount(int scale) {
     int depth = scale - getRootScale();
     if (depth < 0) {
@@ -221,8 +275,9 @@ template <int D> void MWTree<D>::decrementNodeCount(int scale) {
     }
 }
 
-/** @returns Total number of nodes in the tree, at given depth
- * @param[in] depth: Tree depth to count, negative means count _all_ nodes
+/** @returns Total number of nodes in the tree, at given depth (not in use)
+ *
+ * @param[in] depth: Tree depth (0 depth is the coarsest scale) to count.
  */
 template <int D> int MWTree<D>::getNNodesAtDepth(int depth) const {
     int N = 0;
@@ -240,12 +295,13 @@ template <int D> int MWTree<D>::getSizeNodes() const {
     return sizeof(double) * nCoefs / 1024;
 }
 
-/** Find and return the node with the given NodeIndex, const version.
+/** @brief Finds and returns the node pointer with the given \ref NodeIndex, const version.
  *
- * Recursive routine to find and return the node with a given NodeIndex.
- * This routine returns the appropriate Node, or a NULL pointer if
- * the node does not exist, or if it is a GenNode. Recursion starts at the
- * appropriate rootNode. */
+ * @details Recursive routine to find and return the node with a given
+ * NodeIndex.  This routine returns the appropriate Node, or a NULL
+ * pointer if the node does not exist, or if it is a
+ * GenNode. Recursion starts at the appropriate rootNode.
+ */
 template <int D> const MWNode<D> *MWTree<D>::findNode(NodeIndex<D> idx) const {
     if (getRootBox().isPeriodic()) { periodic::index_manipulation<D>(idx, getRootBox().getPeriodic()); }
     int rIdx = getRootBox().getBoxIndex(idx);
@@ -255,12 +311,13 @@ template <int D> const MWNode<D> *MWTree<D>::findNode(NodeIndex<D> idx) const {
     return root.retrieveNodeNoGen(idx);
 }
 
-/** Find and return the node with the given NodeIndex.
+/** @brief Finds and returns the node pointer with the given \ref NodeIndex.
  *
- * Recursive routine to find and return the node with a given NodeIndex.
- * This routine returns the appropriate Node, or a NULL pointer if
- * the node does not exist, or if it is a GenNode. Recursion starts at the
- * appropriate rootNode. */
+ * @details Recursive routine to find and return the node with a given
+ * NodeIndex.  This routine returns the appropriate Node, or a NULL
+ * pointer if the node does not exist, or if it is a
+ * GenNode. Recursion starts at the appropriate rootNode.
+ */
 template <int D> MWNode<D> *MWTree<D>::findNode(NodeIndex<D> idx) {
     if (getRootBox().isPeriodic()) { periodic::index_manipulation<D>(idx, getRootBox().getPeriodic()); }
     int rIdx = getRootBox().getBoxIndex(idx);
@@ -270,11 +327,13 @@ template <int D> MWNode<D> *MWTree<D>::findNode(NodeIndex<D> idx) {
     return root.retrieveNodeNoGen(idx);
 }
 
-/** Find and return the node with the given NodeIndex.
+/** @brief Finds and returns the node reference with the given NodeIndex.
  *
- * This routine ALWAYS returns the node you ask for, and will generate nodes
- * that does not exist. Recursion starts at the appropriate rootNode and
- * decends from this.*/
+ * @details This routine ALWAYS returns the node you ask for. If the
+ * node does not exist, it will be generated by MW
+ * transform. Recursion starts at the appropriate rootNode and descends
+ * from this.
+ */
 template <int D> MWNode<D> &MWTree<D>::getNode(NodeIndex<D> idx) {
     if (getRootBox().isPeriodic()) periodic::index_manipulation<D>(idx, getRootBox().getPeriodic());
 
@@ -289,11 +348,14 @@ template <int D> MWNode<D> &MWTree<D>::getNode(NodeIndex<D> idx) {
     return *out;
 }
 
-/** Find and return the node with the given NodeIndex.
+/** @brief Finds and returns the node with the given NodeIndex.
  *
- * This routine returns the Node you ask for, or the EndNode on
- * the path to the requested node, and will never create or return GenNodes.
- * Recursion starts at the appropriate rootNode and decends from this. */
+ * @details This routine returns the Node you ask for, or the EndNode
+ * on the path to the requested node, if the requested one is deeper
+ * than the leaf node ancestor. It will never create or return
+ * GenNodes.  Recursion starts at the appropriate rootNode and decends
+ * from this.
+ */
 template <int D> MWNode<D> &MWTree<D>::getNodeOrEndNode(NodeIndex<D> idx) {
     if (getRootBox().isPeriodic()) { periodic::index_manipulation<D>(idx, getRootBox().getPeriodic()); }
     MWNode<D> &root = getRootBox().getNode(idx);
@@ -301,11 +363,13 @@ template <int D> MWNode<D> &MWTree<D>::getNodeOrEndNode(NodeIndex<D> idx) {
     return *root.retrieveNodeOrEndNode(idx);
 }
 
-/** Find and return the node with the given NodeIndex.
+/** @brief Finds and returns the node reference with the given NodeIndex. Const version.
  *
- * This routine returns the Node you ask for, or the EndNode on
- * the path to the requested node, and will never create or return GenNodes.
- * Recursion starts at the appropriate rootNode and decends from this. */
+ * @details This routine ALWAYS returns the node you ask for. If the
+ * node does not exist, it will be generated by MW
+ * transform. Recursion starts at the appropriate rootNode and decends
+ * from this.
+ */
 template <int D> const MWNode<D> &MWTree<D>::getNodeOrEndNode(NodeIndex<D> idx) const {
     if (getRootBox().isPeriodic()) { periodic::index_manipulation<D>(idx, getRootBox().getPeriodic()); }
     const MWNode<D> &root = getRootBox().getNode(idx);
@@ -313,11 +377,15 @@ template <int D> const MWNode<D> &MWTree<D>::getNodeOrEndNode(NodeIndex<D> idx)
     return *root.retrieveNodeOrEndNode(idx);
 }
 
-/** Find and return the node at a given depth that contains a given coordinate.
+/** @brief Finds and returns the node at a given depth that contains a given coordinate.
  *
- * This routine ALWAYS returns the node you ask for, and will generate nodes
- * that does not exist. Recursion starts at the appropriate rootNode and
- * decends from this. */
+ * @param[in] depth: requested node depth from root scale.
+ * @param[in] r: coordinates of an arbitrary point in space
+ *
+ * @details This routine ALWAYS returns the node you ask for, and will
+ * generate nodes that do not exist. Recursion starts at the
+ * appropriate rootNode and decends from this.
+ */
 template <int D> MWNode<D> &MWTree<D>::getNode(Coord<D> r, int depth) {
     MWNode<D> &root = getRootBox().getNode(r);
     if (depth >= 0) {
@@ -327,11 +395,15 @@ template <int D> MWNode<D> &MWTree<D>::getNode(Coord<D> r, int depth) {
     }
 }
 
-/** Find and return the node at a given depth that contains a given coordinate.
+/** @brief Finds and returns the node at a given depth that contains a given coordinate.
+ *
+ * @param[in] depth: requested node depth from root scale.
+ * @param[in] r: coordinates of an arbitrary point in space
  *
- * This routine returns the Node you ask for, or the EndNode on
+ * @details This routine returns the Node you ask for, or the EndNode on
  * the path to the requested node, and will never create or return GenNodes.
- * Recursion starts at the appropriate rootNode and decends from this. */
+ * Recursion starts at the appropriate rootNode and decends from this.
+ */
 template <int D> MWNode<D> &MWTree<D>::getNodeOrEndNode(Coord<D> r, int depth) {
 
     if (getRootBox().isPeriodic()) { periodic::coord_manipulation<D>(r, getRootBox().getPeriodic()); }
@@ -340,11 +412,15 @@ template <int D> MWNode<D> &MWTree<D>::getNodeOrEndNode(Coord<D> r, int depth) {
     return *root.retrieveNodeOrEndNode(r, depth);
 }
 
-/** Find and return the node at a given depth that contains a given coordinate.
+/** @brief Finds and returns the node at a given depth that contains a given coordinate. Const version
+ *
+ * @param[in] depth: requested node depth from root scale.
+ * @param[in] r: coordinates of an arbitrary point in space
  *
- * This routine returns the Node you ask for, or the EndNode on
+ * @details This routine returns the Node you ask for, or the EndNode on
  * the path to the requested node, and will never create or return GenNodes.
- * Recursion starts at the appropriate rootNode and decends from this. */
+ * Recursion starts at the appropriate rootNode and decends from this.
+ */
 template <int D> const MWNode<D> &MWTree<D>::getNodeOrEndNode(Coord<D> r, int depth) const {
 
     if (getRootBox().isPeriodic()) { periodic::coord_manipulation<D>(r, getRootBox().getPeriodic()); }
@@ -352,6 +428,11 @@ template <int D> const MWNode<D> &MWTree<D>::getNodeOrEndNode(Coord<D> r, int de
     return *root.retrieveNodeOrEndNode(r, depth);
 }
 
+/** @brief Returns the list of all EndNodes
+ *
+ * @details copies the list of all EndNode pointers into a new vector
+ * and retunrs it.
+ */
 template <int D> MWNodeVector<D> *MWTree<D>::copyEndNodeTable() {
     auto *nVec = new MWNodeVector<D>;
     for (int n = 0; n < getNEndNodes(); n++) {
@@ -361,6 +442,12 @@ template <int D> MWNodeVector<D> *MWTree<D>::copyEndNodeTable() {
     return nVec;
 }
 
+/** @brief Recreate the endNodeTable
+ *
+ * @details the endNodeTable is first deleted and then rebuilt from
+ * scratch. It makes use of the TreeIterator to traverse the tree.
+ * 
+ */
 template <int D> void MWTree<D>::resetEndNodeTable() {
     clearEndNodeTable();
     TreeIterator<D> it(*this, TopDown, Hilbert);
@@ -371,6 +458,7 @@ template <int D> void MWTree<D>::resetEndNodeTable() {
     }
 }
 
+
 template <int D> int MWTree<D>::countBranchNodes(int depth) {
     NOT_IMPLEMENTED_ABORT;
 }
@@ -390,7 +478,7 @@ template <int D> int MWTree<D>::countLeafNodes(int depth) {
     //    return nNodes;
 }
 
-/** Traverse tree and count nodes belonging to this rank. */
+/* Traverse tree and count nodes belonging to this rank. */
 template <int D> int MWTree<D>::countNodes(int depth) {
     NOT_IMPLEMENTED_ABORT;
     //    TreeIterator<D> it(*this);
@@ -407,7 +495,7 @@ template <int D> int MWTree<D>::countNodes(int depth) {
     //    return count;
 }
 
-/** Traverse tree and count nodes with allocated coefficients. */
+/* Traverse tree and count nodes with allocated coefficients. */
 template <int D> int MWTree<D>::countAllocNodes(int depth) {
     NOT_IMPLEMENTED_ABORT;
     //    TreeIterator<D> it(*this);
@@ -424,6 +512,8 @@ template <int D> int MWTree<D>::countAllocNodes(int depth) {
     //    return count;
 }
 
+/** @brief Prints a summary of the tree structure on the output file
+ */
 template <int D> std::ostream &MWTree<D>::print(std::ostream &o) const {
     o << "  square norm: " << this->squareNorm << std::endl;
     o << "  root scale: " << this->getRootScale() << std::endl;
@@ -436,7 +526,11 @@ template <int D> std::ostream &MWTree<D>::print(std::ostream &o) const {
     return o;
 }
 
-/** set values for maxSquareNorm in all nodes  */
+/** @brief sets values for maxSquareNorm in all nodes
+ *
+ * @details it defines the upper bound of the squared norm \f$
+ * ||f||^2_{\ldots} \f$ in this node or its descendents
+ */
 template <int D> void MWTree<D>::makeMaxSquareNorms() {
     NodeBox<D> &rBox = this->getRootBox();
     MWNode<D> **roots = rBox.getNodes();
@@ -446,7 +540,10 @@ template <int D> void MWTree<D>::makeMaxSquareNorms() {
     }
 }
 
-/** gives serialIx of a node from its NodeIndex */
+/** @brief gives serialIx of a node from its NodeIndex
+ *
+ * @details Peter will document this!
+ */
 template <int D> int MWTree<D>::getIx(NodeIndex<D> nIdx) {
     if (this->isLocal == false) MSG_ERROR("getIx only implemented in local representation");
     if(NodeIndex2serialIx.count(nIdx) == 0) return -1;
diff --git a/src/trees/MWTree.h b/src/trees/MWTree.h
index 2bde71767..cd4beb88b 100644
--- a/src/trees/MWTree.h
+++ b/src/trees/MWTree.h
@@ -38,6 +38,27 @@
 
 namespace mrcpp {
 
+/** @class MWTree
+ *
+ * @brief Base class for Multiwavelet tree structures, such as FunctionTree and OperatorTree
+ *
+ * @details The MWTree class is the base class for all tree structures
+ * needed for Multiwavelet calculations. The MWTree is a D-dimensional
+ * tree structure of MWNodes. The tree starts from a set of root nodes
+ * at a common given scale, defining the world box. The most common
+ * settings are either a single root node or \f$ 2^D \f$ root
+ * nodes. Other configurations are however allowed. For example, in 3D
+ * one could have a total of 12 root nodes (a 2x2x3 set of root
+ * nodes). Once the tree structure is generated, each node will have a
+ * parent node (except for the root nodes) and \f$ 2^D \f$ child nodes
+ * (except for leaf nodes). Most of the methods deal with traversing
+ * the tree structure in different ways to fetch specific nodes. Some
+ * of them will return a node present in the tree; some other methods
+ * will generate the required node on the fly using the MW transform;
+ * some methods will return an empty pointer if the node is not
+ * present. See specific methods for details.
+ *
+ */
 template <int D> class MWTree {
 public:
     MWTree(const MultiResolutionAnalysis<D> &mra, const std::string &n);
@@ -58,6 +79,7 @@ template <int D> class MWTree {
     int getKp1_d() const { return this->kp1_d; }
     int getDim() const { return D; }
     int getTDim() const { return (1 << D); }
+    /** @returns the total number of nodes in the tree */
     int getNNodes() const { return getNodeAllocator().getNNodes(); }
     int getNNegScales() const { return this->nodesAtNegativeDepth.size(); }
     int getRootScale() const { return this->rootBox.getScale(); }
@@ -65,6 +87,7 @@ template <int D> class MWTree {
     int getNNodesAtDepth(int i) const;
     int getSizeNodes() const;
 
+    /** @returns */
     NodeBox<D> &getRootBox() { return this->rootBox; }
     const NodeBox<D> &getRootBox() const { return this->rootBox; }
     const MultiResolutionAnalysis<D> &getMRA() const { return this->MRA; }