diff --git a/.clang-tidy b/.clang-tidy
new file mode 100644
index 0000000..a20e2b8
--- /dev/null
+++ b/.clang-tidy
@@ -0,0 +1,20 @@
+Checks: >
+  clang-diagnostic-*,
+  clang-analyzer-*,
+  modernize-*,
+  readability-*,
+  performance-*,
+  bugprone-*,
+  -modernize-use-trailing-return-type,
+  -readability-magic-numbers,
+  -readability-named-parameter,
+  -readability-const-return-type,
+  -readability-identifier-length,
+  -readability-braces-around-statements,
+  -readability-implicit-bool-conversion,
+  -readability-function-cognitive-complexity,
+  -modernize-avoid-c-arrays,
+  -bugprone-easily-swappable-parameters,
+  -modernize-use-equals-default,
+  -readability-convert-member-functions-to-static
+HeaderFilterRegex: '.*(src|test|demo)/.*\.(h|cc)$'
diff --git a/.editorconfig b/.editorconfig
index 4fb9c72..3502e1c 100644
--- a/.editorconfig
+++ b/.editorconfig
@@ -6,7 +6,11 @@ end_of_line = lf
 insert_final_newline = true
 trim_trailing_whitespace = true
 
-[*.{cc,h}]
+[*.{h,cc}]
+indent_style = space
+indent_size = 2
+
+[CMakeLists.txt]
 indent_style = space
 indent_size = 2
 
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 8103405..f55bbaf 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -10,14 +10,8 @@ jobs:
     steps:
       - uses: actions/checkout@v4
 
-      - name: Install GSL
+      - name: Install dependencies
         run: sudo apt-get install -y libgsl-dev
 
-      - name: Configure
-        run: cmake -B build -DCMAKE_BUILD_TYPE=Debug
-
-      - name: Build
-        run: cmake --build build
-
-      - name: Test
-        run: build/dvrlib_main
+      - name: Build and test
+        run: make test
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
new file mode 100644
index 0000000..42e0444
--- /dev/null
+++ b/.github/workflows/docs.yml
@@ -0,0 +1,34 @@
+name: Docs
+
+on:
+  push:
+    branches: [master]
+
+permissions:
+  contents: read
+  pages: write
+  id-token: write
+
+jobs:
+  build-and-deploy:
+    runs-on: ubuntu-latest
+    environment:
+      name: github-pages
+      url: ${{ steps.deploy.outputs.page_url }}
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Install dependencies
+        run: sudo apt-get install -y libgsl-dev doxygen graphviz
+
+      - name: Build docs
+        run: make doc
+
+      - name: Upload Pages artifact
+        uses: actions/upload-pages-artifact@v3
+        with:
+          path: doc/html
+
+      - name: Deploy to GitHub Pages
+        id: deploy
+        uses: actions/deploy-pages@v4
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7be0669..905800d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,18 +25,18 @@ target_link_libraries(dvrlib PUBLIC GSL::gsl GSL::gslcblas)
 target_compile_options(dvrlib PRIVATE -Wall)
 
 # Tests
-add_executable(dvrlib_main
+add_executable(dvrlib_test
     test/main.cc
     test/gsl_wrapper_tests.cc
     test/recon_tests.cc
     test/vdi2048_test.cc
 )
-target_link_libraries(dvrlib_main PRIVATE dvrlib Catch2::Catch2)
-target_compile_options(dvrlib_main PRIVATE -Wall)
+target_link_libraries(dvrlib_test PRIVATE dvrlib Catch2::Catch2)
+target_compile_options(dvrlib_test PRIVATE -Wall)
 
 include(CTest)
 include(Catch)
-catch_discover_tests(dvrlib_main)
+catch_discover_tests(dvrlib_test)
 
 # Demo
 add_executable(dvrlib_demo demo/vdi2048_demo.cc demo/vdi2048.cc)
diff --git a/Doxyfile b/Doxyfile
index ec47687..20f95d9 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -854,7 +854,7 @@ WARN_LOGFILE           =
 # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
 # Note: If this tag is empty the current directory is searched.
 
-INPUT                  =
+INPUT                  = README.md src
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -1007,7 +1007,7 @@ FILTER_SOURCE_PATTERNS =
 # (index.html). This can be useful if you have a project on for instance GitHub
 # and want to reuse the introduction page also for the doxygen output.
 
-USE_MDFILE_AS_MAINPAGE =
+USE_MDFILE_AS_MAINPAGE = README.md
 
 #---------------------------------------------------------------------------
 # Configuration options related to source browsing
diff --git a/Makefile b/Makefile
index 178b02a..961461e 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,8 @@
-.PHONY: all demo test coverage coverage-html clean
+SHELL := /bin/bash
+SRCS = $(shell find src test demo -name '*.h' -o -name '*.cc')
+
+.PHONY: all demo test coverage coverage-html check doc clean
+
 
 BUILD     = build
 BUILD_COV = build-cov
@@ -7,28 +11,35 @@ COV_OUT   = coverage-html
 all: $(BUILD)/CMakeCache.txt
 	cmake --build $(BUILD) -- -s
 
-demo: all
+run-demo: all
 	$(BUILD)/dvrlib_demo
 
-test: all
-	@$(BUILD)/dvrlib_main
+run-demo-zalg: all
+	$(BUILD)/dvrlib_demo --keep-free
+
+run-demo-compare: all
+	@diff --width=160 -y <($(BUILD)/dvrlib_demo) <($(BUILD)/dvrlib_demo --keep-free) || true
+
+run-tests: all
+	@$(BUILD)/dvrlib_test
 
-coverage: $(BUILD_COV)/CMakeCache.txt
+run-coverage: $(BUILD_COV)/CMakeCache.txt
 	cmake --build $(BUILD_COV) -- -s
-	$(BUILD_COV)/dvrlib_main
+	$(BUILD_COV)/dvrlib_test
 	@cd $(BUILD_COV)/CMakeFiles/dvrlib.dir/src && \
 	    gcov *.gcno 2>/dev/null \
 	        | awk '/File.*\/src\/[^/]+\.cc/{f=1; print} f && /Lines executed/{print; f=0}'; \
 	    rm -f *.gcov 2>/dev/null || true
 
-coverage-html: $(BUILD_COV)/CMakeCache.txt
+run-coverage-html: $(BUILD_COV)/CMakeCache.txt
 	@which lcov > /dev/null 2>&1 || { echo "lcov not found. Install with: sudo apt install lcov"; exit 1; }
 	cmake --build $(BUILD_COV) -- -s
-	$(BUILD_COV)/dvrlib_main
+	$(BUILD_COV)/dvrlib_test
 	lcov --capture --directory $(BUILD_COV) --output-file $(BUILD_COV)/coverage.info
 	lcov --remove $(BUILD_COV)/coverage.info '*/build-cov/_deps/*' '/usr/*' --output-file $(BUILD_COV)/coverage.info
 	genhtml $(BUILD_COV)/coverage.info --output-directory $(COV_OUT)
-	@echo "Report: $(COV_OUT)/index.html"
+	@echo "================================================="
+	@echo "Report written to: $(COV_OUT)/index.html"
 
 $(BUILD)/CMakeCache.txt:
 	cmake -B $(BUILD) -S . -DCMAKE_BUILD_TYPE=Debug
@@ -38,5 +49,19 @@ $(BUILD_COV)/CMakeCache.txt:
 	    -DCMAKE_CXX_FLAGS="--coverage" \
 	    -DCMAKE_EXE_LINKER_FLAGS="--coverage"
 
+check: $(BUILD)/CMakeCache.txt
+	# Commented out, clang-tidy reports too much (sometimes dangerous) nonsense
+	#	@echo "--- clang-tidy ---"
+	#	@find src test demo -name '*.cc' | xargs clang-tidy -p $(BUILD)
+	@echo "--- clang-format ---"
+	@clang-format --dry-run --Werror $(SRCS) && echo "format OK"
+
+format_source:
+	@echo "--- clang-format ---"
+	@clang-format $(SRCS) -i --verbose
+
+doc: $(BUILD)/CMakeCache.txt
+	cmake --build $(BUILD) --target doc
+
 clean:
 	rm -rf $(BUILD) $(BUILD_COV) $(COV_OUT)
diff --git a/README.md b/README.md
index dccab91..1d84568 100644
--- a/README.md
+++ b/README.md
@@ -25,24 +25,33 @@ On Debian/Ubuntu:
 sudo apt install libgsl-dev cmake build-essential
 ```
 
-## Building
+## Building and Testing
 
 ```bash
-cmake -B build -DCMAKE_BUILD_TYPE=Debug   # configure
-cmake -B build                            # build library + test binary
-cmake -B build --target doc               # generate API docs (requires Doxygen)
-rm -rf build                              # clean (remove build directory)
+make                   # configure (if needed) and build
+make run-tests         # build and run tests
+make run-demo          # build and run the VDI 2048 demo
+make run-demo-zalg     # same, using the Z-algorithm (keeps free variables)
+make run-demo-compare  # show both variants side by side
+make run-coverage      # build with coverage flags and show per-file line coverage
+make run-coverage-html # same, but generate HTML report (requires lcov)
+make check             # run clang-format check
+make clean             # remove all build directories
 ```
 
-Produces `build/libdvrlib.a` (static library) and `build/dvrlib_main` (test binary).
-
-## Running Tests
+For the HTML coverage report:
 
 ```bash
-cmake -B build && build/dvrlib_main
+sudo apt install lcov
+make run-coverage-html
+# open coverage-html/index.html
 ```
 
-The test binary runs the GSL wrapper tests, the reconciliation tests, and the VDI 2048 example. A clean exit (return code 0) means all tests pass; an assertion failure or uncaught exception indicates a failure.
+To build the API docs (requires Doxygen):
+
+```bash
+make doc
+```
 
 ## Usage Example
 
@@ -80,13 +89,18 @@ int main() {
 
 ```
 src/
-  gsl_wrapper.{h,cc}       GSL vector/matrix RAII wrappers
-  gsl_wrapper_tests.{h,cc} Tests for the wrapper layer
-  recon.{h,cc}              Reconciliation algorithms
-  recon_system.{h,cc}       High-level system description
-  recon_tests.{h,cc}        Tests for reconciliation
-  vdi2048.{h,cc}            VDI 2048 worked example
-  main.cc                   Test runner entry point
+  gsl_wrapper.{h,cc}    GSL vector/matrix RAII wrappers
+  recon.{h,cc}          Reconciliation algorithms
+  recon_system.{h,cc}   High-level system description
+test/
+  main.cc               Catch2 test runner entry point
+  gsl_wrapper_tests.cc  Tests for the wrapper layer
+  recon_tests.cc        Tests for reconciliation
+  vdi2048_test.cc       VDI 2048 regression test
+  custom_reporter.h     Custom Catch2 summary listener
+demo/
+  vdi2048.{h,cc}        VDI 2048 worked example
+  vdi2048_demo.cc       Demo entry point
 ```
 
 ## License
diff --git a/demo/utils.h b/demo/utils.h
index 8f055ba..b55ccf5 100644
--- a/demo/utils.h
+++ b/demo/utils.h
@@ -3,12 +3,21 @@
 
 #include <cstdio>
 #include <sstream>
+#include <string>
 
-#define PRINT(x) \
-  do { std::ostringstream _s; _s << (x); printf(#x ": %s\n", _s.str().c_str()); } while(0)
+#define PRINT(x)                           \
+  do {                                     \
+    std::ostringstream _s;                 \
+    _s << (x);                             \
+    printf(#x ": %s\n", _s.str().c_str()); \
+  } while(0)
 #define PRINT_SEP \
   printf("======================================================================\n")
-#define PRINT_TITLE(x) \
-  do { PRINT_SEP; printf("     %s\n", (x)); PRINT_SEP; } while(0)
+#define PRINT_TITLE(x)                        \
+  do {                                        \
+    PRINT_SEP;                                \
+    printf("     %s\n", (std::string() + x).c_str()); \
+    PRINT_SEP;                                \
+  } while(0)
 
 #endif  // DVRLIB_UTILS_H
diff --git a/demo/vdi2048.cc b/demo/vdi2048.cc
index 8fcb340..875a814 100644
--- a/demo/vdi2048.cc
+++ b/demo/vdi2048.cc
@@ -1,5 +1,4 @@
 #include <iostream>
-#include <vector>
 
 #include "gsl_wrapper.h"
 #include "recon.h"
@@ -9,9 +8,21 @@
 
 using namespace dvrlib;
 
-void example_VDI2048_1(bool increase_var = false, bool keep_free = true) {
+void example_VDI2048(bool keep_free) {
+  PRINT_TITLE("VDI 2048 demo - " + (keep_free ? "z-algorithm" : "original") + " version.");
+  printf("Performing validation and reconsiliation algorithm according\n");
+  printf("to VDI 2048 given the Appendix A.\n\n");
+  if(keep_free) {
+    printf("Using Z-algorithm (i.e. not removing free variables first).\n");
+    printf("Some intermediate quantities are not computed, but aren't necessary\n");
+    printf("for the final result.\n");
+  } else {
+    printf("Removing free variables first, very close to original VDI 2048.\n");
+  }
+
   recon_system system;
 
+  // Define the measured variables (A6) (or A7, just for the values)
   system.add_var("m_FDKeI", 46.241, 0.800);
   system.add_var("m_FDKeII", 45.668, 0.790);
   system.add_var("m_SpI", 44.575, 0.535);
@@ -31,66 +42,74 @@ void example_VDI2048_1(bool increase_var = false, bool keep_free = true) {
     system.add_var("m_HDAnz", 0, -1);
   }
 
+  // Define correlation coefficients (A8)
   system.add_covariance_coeff("m_FDKeI", "m_FDKeII", 0.2);
   system.add_covariance_coeff("m_SpI", "m_SpII", 0.4);
 
-  if(increase_var) {
-    system.change_var("m_FDKeI", 46.241, 2.500);
-    system.change_var("m_FDKeII", 45.668, 2.500);
-  }
+  PRINT_TITLE("Original variables and confidences (A6)");
+  system.print_vars();
 
   // get and print the covariance matrix
   matrix S_x = system.get_covariance_matrix();
-  PRINT_TITLE("Matrix S_X (compare A9, values from A6 and A8)");
+  PRINT_TITLE("Matrix S_X (A9)");
   PRINT(S_x);
 
-  // define the constraints
+  // Define the constraints (A1) - (A4)
   //   m_FD1 = m_FDKeI + m_FDKeII - 0.2 * m_V
   //   m_FD2 = m_SpI + m_SpII - 0.6 * m_V
   //   m_FD3 = m_HK + m_A7 + m_A6 + m_A5 + 0.4 * m_V
   //   m_HDAnz = m_A7 + m_A6 + m_A5
-  // plus (if keep_free == true)
+  // If we keep the free vars (keep_free == true), we also need (A12, three lines)
   //   M_FD1 - M_FD2 = 0
   //   M_FD2 - M_FD3 = 0
   //   M_HDAnz - M_HDNK = 0
+  // With contraints this is just putting the constraints as they are into the matrix
   auto F_with_free = []() -> matrix {
     double Fc[][15] = {
-      {1, 1, 0, 0, -0.2, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0},
-      {0, 0, 1, 1, -0.6, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0},
-      {0, 0, 0, 0, 0.4, 1, 1, 1, 1, 0, 0, 0, 0, -1, 0},
-      {0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, -1},
-      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0},
-      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0},
-      {0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1},
+      {1, 1, 0, 0, -0.2, 0, 0, 0, 0,  0, 0, -1,  0,  0,  0},
+      {0, 0, 1, 1, -0.6, 0, 0, 0, 0,  0, 0,  0, -1,  0,  0},
+      {0, 0, 0, 0,  0.4, 1, 1, 1, 1,  0, 0,  0,  0, -1,  0},
+      {0, 0, 0, 0,  0.0, 0, 1, 1, 1,  0, 0,  0,  0,  0, -1},
+      {0, 0, 0, 0,  0.0, 0, 0, 0, 0,  0, 0,  1, -1,  0,  0},
+      {0, 0, 0, 0,  0.0, 0, 0, 0, 0,  0, 0,  0,  1, -1,  0},
+      {0, 0, 0, 0,  0.0, 0, 0, 0, 0, -1, 0,  0,  0,  0,  1},
     };
-    return matrix(7, 15, Fc);
+    return {7, 15, Fc};
   };
 
+  // Defining the constrained after having solved for the
+  // free variables (A14)
   auto F_without_free = []() -> matrix {
     double Fc[][11] = {
-      {1, 1, -1, -1, 0.4, 0, 0, 0, 0, 0, 0},
-      {0, 0, 1, 1, -1, -1, -1, -1, -1, 0, 0},
-      {0, 0, 0, 0, 0, 0, 1, 1, 1, -1, 0},
+      {1, 1, -1, -1, 0.4,  0,  0,  0,  0,  0, 0},
+      {0, 0,  1,  1,  -1, -1, -1, -1, -1,  0, 0},
+      {0, 0,  0,  0,   0,  0,  1,  1,  1, -1, 0},
     };
-    return matrix(3, 11, Fc);
+    return {3, 11, Fc};
   };
 
   matrix F = keep_free ? F_with_free() : F_without_free();
+  vector x = system.get_values();
 
-  PRINT_TITLE("Constraints (compare A1-A4 and A12)");
-  system.print_constraints(F);
+  if(keep_free) {
+    PRINT_TITLE("Constraints (compare A1-A4 and A12)");
+    system.print_constraints(F);
+  } else {
+    PRINT_TITLE("Constraints (A13 middle)");
+    system.print_constraints(F);
 
-  PRINT_TITLE("Constraints as Matrix F + (+extra cols)");
-  PRINT(F);
+    PRINT_TITLE("Contradictions (A13 right)");
+    PRINT(F * x);
 
-  vector x = system.get_values();
+    PRINT_TITLE("Constraints as Matrix F (A14)");
+    PRINT(F);
 
-  // A15
-  // PRINT_TITLE("Matrix F*S_X (A15)");
-  // PRINT(F*S_x);
+    PRINT_TITLE("Matrix F*S_X (A15)");
+    PRINT(F * S_x);
 
-  PRINT_TITLE("Matrix F*X (A17)");
-  PRINT(F * x);
+    PRINT_TITLE("Matrix F*S_X*F^T (A16)");
+    PRINT(F * S_x * F.transpose());
+  }
 
   PRINT_SEP;
 
@@ -112,6 +131,30 @@ void example_VDI2048_1(bool increase_var = false, bool keep_free = true) {
   PRINT_TITLE("Matrix S_v (compare A21)");
   PRINT(S_v);
 
+  // Changing the confidences of m_FDKeI and m_FDKeII
+  //===============================================================================
+  system.change_var("m_FDKeI", 46.241, 2.500);
+  system.change_var("m_FDKeII", 45.668, 2.500);
+
+  PRINT_TITLE("Original variables and confidences (A24)");
+  system.print_vars();
+
+  S_x = system.get_covariance_matrix();
+  PRINT_TITLE("Matrix S_X (A25)");
+  PRINT(S_x);
+
+  lin_recon(F * x, S_x, F, v);
+
+  PRINT_TITLE("Vector v (compare A30)");
+  PRINT(v);
+
+  PRINT_TITLE("Updated vector x + v (A31)");
+  PRINT(x + v);
+
+  lin_cov_update(S_x, F, S_v);
+  PRINT_TITLE("Matrix S_v (compare A32)");
+  PRINT(S_v);
+
   // compute and print covariance update S_X_new
   matrix S_x_new(S_x.size1(), S_x.size2());
   S_x_new = S_x - S_v;
@@ -128,18 +171,10 @@ void example_VDI2048_1(bool increase_var = false, bool keep_free = true) {
   vector confint_new(S_x.size1());
   extract_confidence(S_x_new, confint_new);
 
-  PRINT_TITLE("confint_new (compare A34)");
-  PRINT(confint_new);
-
-  PRINT_TITLE("Original variables and confidences");
-  system.print_vars();
-
-  PRINT_TITLE("Updated variables and confidences");
+  PRINT_TITLE("Updated variables and confidences (A34)");
   auto updated_sys = system.updated(x_new, confint_new);
   updated_sys.print_vars();
 
+  PRINT_TITLE("New confidences (A34)");
+  PRINT(confint_new);
 }
-
-void example_VDI2048() {
-  example_VDI2048_1(true, true);
-}
\ No newline at end of file
diff --git a/demo/vdi2048.h b/demo/vdi2048.h
index cecf08b..0fd007e 100644
--- a/demo/vdi2048.h
+++ b/demo/vdi2048.h
@@ -1,6 +1,6 @@
 #ifndef DVRLIB_VDI2048_H
 #define DVRLIB_VDI2048_H
 
-void example_VDI2048();
+void example_VDI2048(bool keep_free = true);
 
 #endif  // DVRLIB_VDI2048_H
diff --git a/demo/vdi2048_demo.cc b/demo/vdi2048_demo.cc
index 42a63ba..7078c1c 100644
--- a/demo/vdi2048_demo.cc
+++ b/demo/vdi2048_demo.cc
@@ -1,10 +1,34 @@
 #include "gsl_wrapper.h"
 #include "vdi2048.h"
 
+#include <cstdio>
+#include <cstring>
+
 using namespace dvrlib;
 
-int main() {
-    gsl_enable_exceptions();
-    example_VDI2048();
-    return 0;
+int main(int argc, char* argv[]) {
+  bool keep_free = false;
+
+  for(int i = 1; i < argc; i++) {
+    if(strcmp(argv[i], "--keep-free") == 0) {
+      keep_free = true;
+    } else if(strcmp(argv[i], "--help") == 0 || strcmp(argv[i], "-h") == 0) {
+      printf("Usage: %s [options]\n", argv[0]);
+      printf("\n");
+      printf("Options:\n");
+      printf("  --keep-free   Keep free variables in the system of equations\n");
+      printf("                (Z-algorithm). Default: eliminate free variables first.\n");
+      printf("  -h, --help    Show this help message.\n");
+      return 0;
+    } else {
+      printf("Unknown option: %s\n", argv[i]);
+      printf("Run with --help for usage.\n");
+      return 1;
+    }
+  }
+
+  gsl_enable_exceptions();
+  set_vdi_print_format(true);
+  example_VDI2048(keep_free);
+  return 0;
 }
diff --git a/src/dvr_assert.h b/src/dvr_assert.h
index 6e4ed50..f70d39d 100644
--- a/src/dvr_assert.h
+++ b/src/dvr_assert.h
@@ -9,15 +9,16 @@ namespace dvrlib {
 struct assertion_error : public std::logic_error {
   assertion_error(const char* cond, const char* file, int line)
       : std::logic_error(std::string("assertion failed: ") + cond +
-                         " (" + file + ":" + std::to_string(line) + ")") {}
+                         " (" + file + ":" + std::to_string(line) + ")") {
+  }
 };
 
 }  // namespace dvrlib
 
-#define dvr_assert(cond) \
-  do { \
-    if (!(cond)) \
+#define dvr_assert(cond)                                        \
+  do {                                                          \
+    if(!(cond))                                                 \
       throw dvrlib::assertion_error(#cond, __FILE__, __LINE__); \
-  } while (0)
+  } while(0)
 
 #endif  // DVRLIB_DVR_ASSERT_H
diff --git a/src/gsl_wrapper.cc b/src/gsl_wrapper.cc
index ec3ca1f..ff5cc41 100644
--- a/src/gsl_wrapper.cc
+++ b/src/gsl_wrapper.cc
@@ -1,7 +1,7 @@
 #include "gsl_wrapper.h"
 #include "dvr_assert.h"
 
-#include <cmath>
+#include <iomanip>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_linalg.h>
@@ -20,12 +20,41 @@ void gsl_enable_exceptions() {
   gsl_set_error_handler(gsl_err_handler);
 }
 
+static bool g_vdi_format = false;
+
+void set_vdi_print_format(bool enable) {
+  g_vdi_format = enable;
+}
+
+class ostream_format_guard {
+public:
+  ostream_format_guard(std::ostream& out, std::streamsize precision)
+      : out_(out),
+        precision_(out.precision()),
+        flags_(out.flags()),
+        fill_(out.fill()) {
+    out_.precision(precision);
+    out_.setf(std::ios::fixed, std::ios::floatfield);
+    // out_.setf(std::ios::internal, std::ios::adjustfield);
+    out_.fill(' ');
+  }
+
+  ~ostream_format_guard() {
+    out_.precision(precision_);
+    out_.flags(flags_);
+    out_.fill(fill_);
+  }
+
+private:
+  std::ostream& out_;
+  std::streamsize precision_;
+  std::ios::fmtflags flags_;
+  char fill_;
+};
+
 ////////////////////////////////////////////////////////////////////////////
 // vector class
 ////////////////////////////////////////////////////////////////////////////
-vector::vector() {
-  v = 0;  // private ctor only for class vector_view
-}
 
 vector::vector(int n) {
   v = gsl_vector_alloc(n);
@@ -49,8 +78,7 @@ vector::vector(const vector& src) {
 }
 
 vector::~vector() {
-  if(v && v->owner)
-    gsl_vector_free(v);
+  gsl_vector_free(v);
 }
 
 gsl_vector* vector::gsl_internal() {
@@ -62,7 +90,7 @@ const gsl_vector* vector::gsl_internal() const {
 }
 
 int vector::size() const {
-  return v->size;
+  return static_cast<int>(v->size);
 }
 
 void vector::set(int i, double val) {
@@ -84,6 +112,11 @@ vector& vector::operator=(const vector& src) {
   return *this;
 }
 
+vector& vector::operator=(const vector_view& src) {
+  gsl_vector_memcpy(v, src.gsl_internal());
+  return *this;
+}
+
 vector& vector::operator+=(const vector& src) {
   gsl_vector_add(v, src.v);
   return *this;
@@ -145,7 +178,7 @@ double vector::norm2() const {
   return gsl_blas_dnrm2(v);
 }
 
-std::ostream& operator<<(std::ostream& out, const vector& vec) {
+std::ostream& print_std_format(std::ostream& out, const vector& vec) {
   out << "[";
   for(int i = 0; i < vec.size(); i++) {
     if(i)
@@ -156,46 +189,85 @@ std::ostream& operator<<(std::ostream& out, const vector& vec) {
   return out;
 }
 
+std::ostream& print_vdi_format(std::ostream& out, const vector& vec) {
+  ostream_format_guard format(out, 3);
+
+  out << "[";
+  for(int i = 0; i < vec.size(); i++) {
+    if(i)
+      out << ", ";
+    auto value = vec.get(i);
+    if(std::abs(value) < 1e-4)
+      out << "    0  ";
+    else
+      out << std::setw(7) << value;
+  }
+  out << "]";
+  return out;
+}
+
+std::ostream& operator<<(std::ostream& out, const vector& vec) {
+  if(g_vdi_format)
+    return print_vdi_format(out, vec);
+  else
+    return print_std_format(out, vec);
+}
+
 ////////////////////////////////////////////////////////////////////////////
 // vector_view class
 ////////////////////////////////////////////////////////////////////////////
 
-vector_view::vector_view(gsl_vector_view _vv) {
+vector_view::vector_view(const gsl_vector_view& _vv) {
   vv = _vv;
-  v = &vv.vector;
 }
 
 vector_view::vector_view(const vector_view& src) {
   vv = src.vv;
-  v = &vv.vector;
 }
 
-gsl_vector_view* vector_view::gsl_internal() {
-  return &vv;
+int vector_view::size() const {
+  return static_cast<int>(vv.vector.size);
+}
+
+double vector_view::get(int i) const {
+  return gsl_vector_get(&vv.vector, i);
+}
+
+void vector_view::set(int i, double val) {
+  gsl_vector_set(&vv.vector, i, val);
+}
+
+const gsl_vector* vector_view::gsl_internal() const {
+  return &vv.vector;
 }
 
 vector_view& vector_view::operator=(const vector& src) {
-  if(&src == this)
-    return *this;
-  gsl_vector_memcpy(v, src.v);
+  gsl_vector_memcpy(&vv.vector, src.gsl_internal());
+  return *this;
+}
+
+vector_view& vector_view::operator=(const vector_view& src) {
+  gsl_vector_memcpy(&vv.vector, &src.vv.vector);
   return *this;
 }
 
+vector_view::operator vector() const {
+  vector result(size());
+  gsl_vector_memcpy(result.gsl_internal(), &vv.vector);
+  return result;
+}
+
 ////////////////////////////////////////////////////////////////////////////
 // matrix class
 ////////////////////////////////////////////////////////////////////////////
 
-matrix::matrix() {
-  m = 0;  // private constructor for matrix_view only
-}
-
 matrix::matrix(int n1, int n2, bool id, const double* diag) {
   m = gsl_matrix_alloc(n1, n2);
   if(!id)
     gsl_matrix_set_zero(m);
   else
     gsl_matrix_set_identity(m);
-  if(diag) {
+  if(diag != nullptr) {
     for(int i = 0; i < n1 && i < n2; i++)
       set(i, i, diag[i]);
   }
@@ -218,8 +290,7 @@ matrix::matrix(const gsl_matrix* src) {
 }
 
 matrix::~matrix() {
-  if(m && m->owner)
-    gsl_matrix_free(m);
+  gsl_matrix_free(m);
 }
 
 gsl_matrix* matrix::gsl_internal() {
@@ -231,11 +302,11 @@ const gsl_matrix* matrix::gsl_internal() const {
 }
 
 int matrix::size1() const {
-  return m->size1;
+  return static_cast<int>(m->size1);
 }
 
 int matrix::size2() const {
-  return m->size2;
+  return static_cast<int>(m->size2);
 }
 
 void matrix::set(int i, int j, double val) {
@@ -263,6 +334,11 @@ matrix& matrix::operator=(const matrix& src) {
   return *this;
 }
 
+matrix& matrix::operator=(const matrix_view& src) {
+  gsl_matrix_memcpy(m, src.gsl_internal());
+  return *this;
+}
+
 matrix matrix::operator+(const matrix& src) const {
   matrix c(*this);
   dvr_assert(c.size1() == src.size1());
@@ -340,7 +416,8 @@ struct permutation {
 };
 
 matrix matrix::inverse() const {
-  matrix LU(*this), B(*this);
+  matrix LU(*this);
+  matrix B(*this);
   permutation p(size1());
   int sign;
   gsl_linalg_LU_decomp(LU.m, p.p, &sign);
@@ -393,7 +470,7 @@ const matrix_view matrix::submatrix(int k1, int k2, int n1, int n2) const {
   return mv;
 }
 
-std::ostream& operator<<(std::ostream& out, const matrix& mat) {
+std::ostream& print_std_format(std::ostream& out, const matrix& mat) {
   for(int i = 0; i < mat.size1(); i++) {
     if(i)
       out << "]" << std::endl
@@ -410,31 +487,86 @@ std::ostream& operator<<(std::ostream& out, const matrix& mat) {
   return out;
 }
 
+std::ostream& print_vdi_format(std::ostream& out, const matrix& mat) {
+  ostream_format_guard format(out, 4);
+
+  for(int i = 0; i < mat.size1(); i++) {
+    if(!i)
+      out << " [";
+    out << std::endl
+        << "  ";
+
+    for(int j = 0; j < mat.size2(); j++) {
+      if(j)
+        out << " ";
+
+      auto value = mat.get(i, j);
+      if(std::abs(value) < 1e-4)
+        out << "    0  ";
+      else
+        out << std::setw(7) << value;
+    }
+  }
+  out << std::endl
+      << "]";
+  return out;
+}
+
+std::ostream& operator<<(std::ostream& out, const matrix& mat) {
+  if(g_vdi_format)
+    return print_vdi_format(out, mat);
+  else
+    return print_std_format(out, mat);
+}
+
 ////////////////////////////////////////////////////////////////////////////
 // matrix_view class
 ////////////////////////////////////////////////////////////////////////////
 
-matrix_view::matrix_view(gsl_matrix_view _mv) {
+matrix_view::matrix_view(const gsl_matrix_view& _mv) {
   mv = _mv;
-  m = &mv.matrix;
 }
 
 matrix_view::matrix_view(const matrix_view& src) {
   mv = src.mv;
-  m = &mv.matrix;
 }
 
-gsl_matrix_view* matrix_view::gsl_internal() {
-  return &mv;
+int matrix_view::size1() const {
+  return static_cast<int>(mv.matrix.size1);
+}
+
+int matrix_view::size2() const {
+  return static_cast<int>(mv.matrix.size2);
+}
+
+double matrix_view::get(int i, int j) const {
+  return gsl_matrix_get(&mv.matrix, i, j);
+}
+
+void matrix_view::set(int i, int j, double val) {
+  gsl_matrix_set(&mv.matrix, i, j, val);
+}
+
+const gsl_matrix* matrix_view::gsl_internal() const {
+  return &mv.matrix;
 }
 
 matrix_view& matrix_view::operator=(const matrix& src) {
-  if(&src == this)
-    return *this;
-  gsl_matrix_memcpy(m, src.m);
+  gsl_matrix_memcpy(&mv.matrix, src.gsl_internal());
   return *this;
 }
 
+matrix_view& matrix_view::operator=(const matrix_view& src) {
+  gsl_matrix_memcpy(&mv.matrix, &src.mv.matrix);
+  return *this;
+}
+
+matrix_view::operator matrix() const {
+  matrix result(size1(), size2());
+  gsl_matrix_memcpy(result.gsl_internal(), &mv.matrix);
+  return result;
+}
+
 matrix operator*(double d, const matrix& src) {
   return src * d;
 }
@@ -444,9 +576,9 @@ matrix operator*(double d, const matrix& src) {
 ////////////////////////////////////////////////////////////////////////////
 
 #define PRINT_START(cls) \
-  out << cls << " {" << std::endl
+  out << (cls) << " {" << std::endl
 #define PRINT_FIELD(s, field) \
-  out << "  " << #field << ": " << s.field << std::endl
+  out << "  " << #field << ": " << (s).field << std::endl
 #define PRINT_END() \
   out << "}" << std::endl
 
diff --git a/src/gsl_wrapper.h b/src/gsl_wrapper.h
index 14d4fdd..c2af73a 100644
--- a/src/gsl_wrapper.h
+++ b/src/gsl_wrapper.h
@@ -9,65 +9,130 @@
 
 namespace dvrlib {
 
+/** Exception type thrown by GSL when \c gsl_enable_exceptions() is active. */
 struct gsl_exception {
-  const char* reason;
-  const char* file;
-  int line;
-  int gsl_errno;
+  const char* reason;  ///< Human-readable error description
+  const char* file;    ///< Source file where the error occurred
+  int line;            ///< Line number where the error occurred
+  int gsl_errno;       ///< GSL error code (see \c gsl_errno.h)
 };
 
+/** Install a GSL error handler that throws \c gsl_exception instead of aborting. */
 void gsl_enable_exceptions();
 
+/**
+ * Switch between standard and VDI 2048 pretty-print format for vectors and
+ * matrices written to \c std::ostream.
+ *
+ * @param[in] enable \c true selects VDI format, \c false selects standard format.
+ */
+void set_vdi_print_format(bool enable);
+
 class vector_view;
 
+/**
+ * RAII wrapper around a heap-allocated \c gsl_vector.
+ *
+ * Owns the underlying GSL allocation and frees it on destruction.
+ * Arithmetic operators follow GSL semantics (element-wise).
+ */
 class vector {
   gsl_vector* v;
-  vector();
 
 public:
+  /** Allocate an uninitialised vector of length \c n. */
   vector(int n);
+  /** Allocate a vector of length \c n with all elements set to \c x. */
   vector(int n, double x);
+  /** Allocate a vector of length \c n and copy values from the array \c x. */
   vector(int n, const double* x);
+  /** Copy constructor. */
   vector(const vector& src);
+  /** Destructor — frees the underlying GSL allocation. */
   ~vector();
 
+  /** Return a pointer to the underlying \c gsl_vector (mutable). */
   gsl_vector* gsl_internal();
+  /** Return a pointer to the underlying \c gsl_vector (read-only). */
   const gsl_vector* gsl_internal() const;
 
+  /** Return the number of elements. */
   int size() const;
+  /** Set element \c i to \c val. */
   void set(int i, double val);
+  /** Return element \c i. */
   double get(int i) const;
+  /** Return element \c i (read-only subscript). */
   double operator[](int i) const;
 
+  /** Copy-assign from another vector. */
   vector& operator=(const vector& src);
+  /** Copy-assign from a view. */
+  vector& operator=(const vector_view& src);
   vector& operator+=(const vector& src);
   vector operator+(const vector& src) const;
   vector& operator-=(const vector& src);
   vector operator-(const vector& src) const;
+  /** Unary negation — returns \c -(*this). */
   vector operator-() const;
 
+  /** Scale all elements by \c d in place. */
   vector& operator*=(double d);
+  /** Return a copy scaled by \c d. */
   vector operator*(double d) const;
+  /** Return the L1 norm (sum of absolute values). */
   double norm1() const;
+  /** Return the L2 (Euclidean) norm. */
   double norm2() const;
 
+  /**
+   * Return a view of \c n elements starting at index \c k.
+   *
+   * @param[in] k Start index (0-based).
+   * @param[in] n Number of elements in the sub-vector.
+   */
   vector_view subvector(int k, int n);
+  /** @copydoc vector::subvector(int,int) */
   const vector_view subvector(int k, int n) const;
-  friend class vector_view;
   friend class matrix;
 };
 
+/** Scalar multiplication \c d*src (commutative with \c vector::operator*). */
 vector operator*(double d, const vector& src);
+/** Write a vector to an output stream. Format depends on \c set_vdi_print_format(). */
 std::ostream& operator<<(std::ostream& out, const vector& vec);
 
-class vector_view : public vector {
+/**
+ * Non-owning view into a contiguous slice of a \c vector or \c matrix row/column.
+ *
+ * Wraps a \c gsl_vector_view. The referenced data is owned by the original
+ * \c vector or \c matrix — the view must not outlive it.
+ * Assigning to a view modifies the underlying data in place.
+ * Use \c operator vector() to obtain an independent copy.
+ */
+class vector_view {
   gsl_vector_view vv;
-  vector_view(gsl_vector_view vv);
+  vector_view(const gsl_vector_view& vv);
 
 public:
+  /** Copy constructor (shallow — both views reference the same data). */
   vector_view(const vector_view& src);
-  gsl_vector_view* gsl_internal();
+
+  /** Return the number of elements. */
+  int size() const;
+  /** Return element \c i. */
+  double get(int i) const;
+  /** Set element \c i to \c val (modifies the underlying data). */
+  void set(int i, double val);
+  /** Return a read-only pointer to the underlying \c gsl_vector. */
+  const gsl_vector* gsl_internal() const;
+
+  /** Copy values from \c src into the viewed memory. */
   vector_view& operator=(const vector& src);
+  /** Copy values from another view into the viewed memory. */
+  vector_view& operator=(const vector_view& src);
+  /** Create an independent \c vector copy of this view. */
+  operator vector() const;
 
   friend class vector;
   friend class matrix;
@@ -75,50 +140,114 @@ class vector_view : public vector {
 
 class matrix_view;
 
+/**
+ * RAII wrapper around a heap-allocated \c gsl_matrix.
+ *
+ * Owns the underlying GSL allocation and frees it on destruction.
+ * Row access via \c operator[] returns a \c vector_view into the matrix data.
+ */
 class matrix {
   gsl_matrix* m;
-  matrix();
 
 public:
-  matrix(int n1, int n2, bool id = false, const double* diag = 0);
+  /**
+   * Allocate an \c n1 × \c n2 matrix.
+   *
+   * @param[in] n1   Number of rows.
+   * @param[in] n2   Number of columns.
+   * @param[in] id   If \c true, initialise as identity (requires \c n1 == \c n2
+   *                 when \c diag is null).
+   * @param[in] diag Optional array of \c n1 diagonal values; if non-null
+   *                 overrides \c id.
+   */
+  matrix(int n1, int n2, bool id = false, const double* diag = nullptr);
+  /** Allocate an \c n1 × \c n2 matrix and copy values from the row-major array \c x. */
   matrix(int n1, int n2, const double* x);
+  /**
+   * Allocate an \c n1 × \c n2 matrix from a 2-D C array.
+   *
+   * The template parameter \c n must match \c n2.
+   */
   template <int n>
   matrix(int n1, int n2, const double (*x)[n]);
 
+  /** Copy constructor. */
   matrix(const matrix& src);
+  /** Construct from an existing (non-owned) \c gsl_matrix pointer. */
   matrix(const gsl_matrix* src);
+  /** Destructor — frees the underlying GSL allocation. */
   ~matrix();
 
+  /** Return a pointer to the underlying \c gsl_matrix (mutable). */
   gsl_matrix* gsl_internal();
+  /** Return a pointer to the underlying \c gsl_matrix (read-only). */
   const gsl_matrix* gsl_internal() const;
 
+  /** Return the number of rows. */
   int size1() const;
+  /** Return the number of columns. */
   int size2() const;
+  /** Set element (\c i, \c j) to \c val. */
   void set(int i, int j, double val);
+  /** Return element (\c i, \c j). */
   double get(int i, int j) const;
+  /** Return a mutable view of row \c i. */
   vector_view operator[](int i);
+  /** Return a read-only view of row \c i. */
   const vector_view operator[](int i) const;
 
+  /** Copy-assign from another matrix. */
   matrix& operator=(const matrix& src);
+  /** Copy-assign from a view. */
+  matrix& operator=(const matrix_view& src);
 
   matrix operator+(const matrix& src) const;
   matrix& operator+=(const matrix& src);
   matrix operator-(const matrix& src) const;
   matrix& operator-=(const matrix& src);
+  /** Unary negation — returns \c -(*this). */
   matrix operator-() const;
 
+  /** Matrix–vector product. */
   vector operator*(const vector& src) const;
+  /** Matrix–matrix product. */
   matrix operator*(const matrix& src) const;
+  /** Return a copy scaled by \c d. */
   matrix operator*(double d) const;
+  /** Scale all elements by \c d in place. */
   matrix& operator*=(double d);
 
+  /** Return the transpose. */
   matrix transpose() const;
+  /** Return the inverse (via LU decomposition). */
   matrix inverse() const;
+  /**
+   * Solve the linear system \c A*x = b and return \c x.
+   *
+   * @param[in] b Right-hand side vector (length \c size1()).
+   */
   vector linsolve(const vector& b) const;
+  /**
+   * Compute the full singular value decomposition \c A = U * diag(s) * V^T.
+   *
+   * @param[out] U Left singular vectors (same size as \c *this).
+   * @param[out] V Right singular vectors (\c size2() × \c size2()).
+   * @param[out] s Singular values (length \c size2()).
+   */
   void svd(matrix& U, matrix& V, vector& s) const;
+  /** Compute and return only the singular values. */
   vector svd() const;
 
+  /**
+   * Return a view of the \c n1 × \c n2 sub-matrix starting at (\c k1, \c k2).
+   *
+   * @param[in] k1 Starting row index (0-based).
+   * @param[in] k2 Starting column index (0-based).
+   * @param[in] n1 Number of rows in the sub-matrix.
+   * @param[in] n2 Number of columns in the sub-matrix.
+   */
   matrix_view submatrix(int k1, int k2, int n1, int n2);
+  /** @copydoc matrix::submatrix(int,int,int,int) */
   const matrix_view submatrix(int k1, int k2, int n1, int n2) const;
   friend class matrix_view;
 };
@@ -131,18 +260,45 @@ matrix::matrix(int n1, int n2, const double (*x)[n]) {
   gsl_matrix_memcpy(m, &src.matrix);
 }
 
+/** Scalar multiplication \c d*src (commutative with \c matrix::operator*). */
 matrix operator*(double d, const matrix& src);
 
-std::ostream& operator<<(std::ostream& out, const matrix& vec);
-
-class matrix_view : public matrix {
+/** Write a matrix to an output stream. Format depends on \c set_vdi_print_format(). */
+std::ostream& operator<<(std::ostream& out, const matrix& mat);
+
+/**
+ * Non-owning view into a rectangular sub-region of a \c matrix.
+ *
+ * Wraps a \c gsl_matrix_view. The referenced data is owned by the original
+ * \c matrix — the view must not outlive it.
+ * Assigning to a view modifies the underlying data in place.
+ * Use \c operator matrix() to obtain an independent copy.
+ */
+class matrix_view {
   gsl_matrix_view mv;
-  matrix_view(gsl_matrix_view mv);
+  matrix_view(const gsl_matrix_view& mv);
 
 public:
+  /** Copy constructor (shallow — both views reference the same data). */
   matrix_view(const matrix_view& src);
-  gsl_matrix_view* gsl_internal();
+
+  /** Return the number of rows. */
+  int size1() const;
+  /** Return the number of columns. */
+  int size2() const;
+  /** Return element (\c i, \c j). */
+  double get(int i, int j) const;
+  /** Set element (\c i, \c j) to \c val (modifies the underlying data). */
+  void set(int i, int j, double val);
+  /** Return a read-only pointer to the underlying \c gsl_matrix. */
+  const gsl_matrix* gsl_internal() const;
+
+  /** Copy values from \c src into the viewed memory. */
   matrix_view& operator=(const matrix& src);
+  /** Copy values from another view into the viewed memory. */
+  matrix_view& operator=(const matrix_view& src);
+  /** Create an independent \c matrix copy of this view. */
+  operator matrix() const;
 
   friend class matrix;
 };
diff --git a/src/recon.cc b/src/recon.cc
index 9bfe549..fbfe0c5 100644
--- a/src/recon.cc
+++ b/src/recon.cc
@@ -161,12 +161,14 @@ int recon(const vector& x,
   return 0;
 }
 
+static constexpr auto z_alpha_95 = 1.959963984540054;
+
 double confint2var(double confint) {
-  return pow((confint / 1.96), 2);
+  return pow((confint / z_alpha_95), 2);
 }
 
 double var2confint(double var) {
-  return 1.96 * sqrt(var);
+  return z_alpha_95 * sqrt(var);
 }
 
 void extract_confidence(const matrix& S_xnew,
diff --git a/src/recon.h b/src/recon.h
index 996e728..004a493 100644
--- a/src/recon.h
+++ b/src/recon.h
@@ -6,27 +6,85 @@ namespace dvrlib {
 class vector;
 class matrix;
 
+/**
+ * Abstract function object mapping \c argtype to \c restype.
+ *
+ * Derive from this to supply residual functions and Jacobians to \c recon().
+ */
 template <class argtype, class restype>
 class func {
 public:
+  virtual ~func() = default;
+  /** Evaluate the function at \c arg. */
   virtual restype operator()(const argtype& arg) const = 0;
 };
 
+/**
+ * Propagate measurement covariance through a linear constraint matrix.
+ *
+ * Computes \c S_v = F * S_x * F^T, the covariance of the reconciled
+ * values given the measurement covariance \c S_x and the constraint
+ * matrix \c F.
+ *
+ * @param[in]  S_x Measurement covariance matrix.
+ * @param[in]  F   Linear constraint (Jacobian) matrix.
+ * @param[out] S_v Resulting covariance of the reconciled values.
+ */
 void lin_cov_update(const matrix& S_x,
                     const matrix& F,
                     matrix& S_v);
 
+/**
+ * Solve the linear reconciliation problem.
+ *
+ * Computes the reconciled values \c v that satisfy the linear constraints
+ * encoded in \c F, given raw measurements \c r and measurement covariance
+ * \c S_x.
+ *
+ * @param[in]  r   Raw measurement vector.
+ * @param[in]  S_x Measurement covariance matrix.
+ * @param[in]  F   Linear constraint matrix.
+ * @param[out] v   Reconciled values.
+ */
 void lin_recon(const vector& r,
                const matrix& S_x,
                const matrix& F,
                vector& v);
 
+/**
+ * Compute the reconciliation correction step for non-linear iteration.
+ *
+ * Given the current iterate \c v, computes the Newton-like correction
+ * \c dv to reduce the residual of the constraint equations.
+ *
+ * @param[in]  r       Raw measurement vector.
+ * @param[in]  S_x_inv Inverse of the measurement covariance matrix.
+ * @param[in]  F       Constraint Jacobian evaluated at the current iterate.
+ * @param[in]  v       Current reconciled values.
+ * @param[out] dv      Correction to apply: next iterate is \c v + dv.
+ */
 void lin_recon_update(const vector& r,
                       const matrix& S_x_inv,
                       const matrix& F,
                       const vector& v,
                       vector& dv);
 
+/**
+ * Solve the non-linear reconciliation problem by iterative linearisation.
+ *
+ * Iterates a Newton-type update until the correction \c dv is smaller than
+ * \c eps (in the L2 norm) or \c maxiter iterations are reached.
+ *
+ * @param[in]  x       Raw measurement vector.
+ * @param[in]  S_x     Measurement covariance matrix.
+ * @param[in]  f       Constraint residual function \c f(v) = 0.
+ * @param[in]  J       Jacobian of \c f with respect to \c v.
+ * @param[out] v       Reconciled values.
+ * @param[out] S_v     Covariance of the reconciled values.
+ * @param[in]  eps     Convergence tolerance on the correction norm (default 1e-6).
+ * @param[in]  maxiter Maximum number of iterations (default 50).
+ * @return Number of iterations performed, or -1 if the method did not converge.
+ */
 int recon(const vector& x,
           const matrix& S_x,
           const func<vector, vector>& f,
@@ -36,17 +94,29 @@ int recon(const vector& x,
           double eps = 1e-6,
           int maxiter = 50);
 
+/**
+ * Extract 95% confidence intervals from a covariance matrix diagonal.
+ *
+ * @param[in]  S_xnew       Covariance matrix of the reconciled values.
+ * @param[out] conf_results Confidence intervals (one per variable).
+ */
 void extract_confidence(const matrix& S_xnew,
                         vector& conf_results);
 
 /**
-   Convert a 95% confidence interval into a variance.
-*/
+ * Convert a 95% confidence interval into a variance.
+ *
+ * Uses the relation \c confint = z_{0.975} * sqrt(var) where
+ * \c z_{0.975} = Φ⁻¹(0.975) ≈ 1.96.
+ */
 double confint2var(double confint);
 
 /**
-   Convert a variance into a 95% confidence interval.
-*/
+ * Convert a variance into a 95% confidence interval.
+ *
+ * Uses the relation \c confint = z_{0.975} * sqrt(var) where
+ * \c z_{0.975} = Φ⁻¹(0.975) ≈ 1.96.
+ */
 double var2confint(double var);
 
 }  // namespace dvrlib
diff --git a/src/recon_system.cc b/src/recon_system.cc
index a37eca7..e994d24 100644
--- a/src/recon_system.cc
+++ b/src/recon_system.cc
@@ -23,7 +23,7 @@ void recon_system::add_covariance_coeff(const char* name1, const char* name2,
 }
 
 int recon_system::find_var(const string& str) const {
-  int n = vars.size();
+  int n = static_cast<int>(vars.size());
   for(int i = 0; i < n; i++) {
     if(vars[i].name == str)
       return i;
@@ -39,9 +39,9 @@ void recon_system::change_var(const char* name, double val, double confint) {
 }
 
 int recon_system::get_number_measured() const {
-  int count = 0, n = vars.size();
-  for(int i = 0; i < n; i++) {
-    if(vars[i].confint >= 0)
+  int count = 0;
+  for(const auto& var : vars) {
+    if(var.confint >= 0)
       count++;
   }
   return count;
@@ -58,14 +58,14 @@ matrix recon_system::get_covariance_matrix() const {
     S_x.set(i, i, confint2var(vars[i].confint));
   }
 
-  for(unsigned int k = 0; k < extra_covs.size(); k++) {
-    int i = find_var(extra_covs[k].var1);
-    int j = find_var(extra_covs[k].var2);
+  for(const auto& extra_cov : extra_covs) {
+    int i = find_var(extra_cov.var1);
+    int j = find_var(extra_cov.var2);
     dvr_assert(i >= 0);
     dvr_assert(j >= 0);
     dvr_assert(i < n);
     dvr_assert(j < n);
-    double rho = extra_covs[k].cov_coeff;
+    double rho = extra_cov.cov_coeff;
     double cov_ii = S_x.get(i, i);
     double cov_jj = S_x.get(j, j);
     double cov_ij = rho * sqrt(cov_ii * cov_jj);
@@ -76,7 +76,7 @@ matrix recon_system::get_covariance_matrix() const {
 }
 
 vector recon_system::get_values() const {
-  int n = vars.size();
+  int n = static_cast<int>(vars.size());
   vector x(n);
   for(int i = 0; i < n; i++) {
     x.set(i, vars[i].value);
@@ -86,13 +86,13 @@ vector recon_system::get_values() const {
 
 recon_system recon_system::updated(const vector& values,
                                    const vector& confints) const {
-  dvr_assert(values.size() == (int)vars.size());
+  dvr_assert(values.size() == static_cast<int>(vars.size()));
   dvr_assert(confints.size() == get_number_measured());
   recon_system result(*this);
-  for(int i = 0; i < (int)vars.size(); i++) {
+  for(int i = 0; i < static_cast<int>(vars.size()); i++) {
     result.vars[i].value = values.get(i);
   }
-  for(int i = 0; i < (int)confints.size(); i++) {
+  for(int i = 0; i < confints.size(); i++) {
     result.vars[i].confint = confints.get(i);
   }
   return result;
@@ -100,12 +100,12 @@ recon_system recon_system::updated(const vector& values,
 
 void recon_system::print_vars() const {
   printf("%-12s%10s%10s\n", "name", "value", "confint");
-  for(unsigned int i = 0; i < vars.size(); i++) {
-    printf("%-12s%10.3f", vars[i].name.c_str(), vars[i].value);
-    if(vars[i].confint < 0)
+  for(const auto& var : vars) {
+    printf("%-12s%10.3f", var.name.c_str(), var.value);
+    if(var.confint < 0)
       printf("    (free)");
     else
-      printf("%10.3f", vars[i].confint);
+      printf("%10.3f", var.confint);
     printf("\n");
   }
 }
diff --git a/src/recon_system.h b/src/recon_system.h
index 16522a7..0fed6df 100644
--- a/src/recon_system.h
+++ b/src/recon_system.h
@@ -8,9 +8,9 @@
 namespace dvrlib {
 
 /**
-   Encapsulates a system to be reconciliated, including variables,
-   confidence intervals, covariance coefficient etc. Constraints still
-   need to be implemented. Also fixed variables should be treated
+   Encapsulates a system to be reconciled, including variables,
+   confidence intervals, covariance coefficient, etc. Constraints still
+   need to be implemented. Also, fixed variables should be treated
    here.
  */
 class recon_system {
@@ -30,8 +30,8 @@ class recon_system {
 
 public:
   /**
-   * Add a variable to be reconciliated to the system. If \c confint is
-   * negative it is assumed that this variable is free (i.e. no
+   * Add a variable to be reconciled to the system. If \c confint is
+   * negative, it is assumed that this variable is free (i.e., no
    * measurements). \todo If \c confint is zero it is assumed that this
    * variable is fixed.
    *
@@ -40,16 +40,69 @@ class recon_system {
    * @param[in] confint The confidence interval
    */
   void add_var(const char* name, double val, double confint);
+
+  /**
+   * Register a covariance coefficient between two variables.
+   *
+   * The coefficient is used when building the covariance matrix via
+   * \c get_covariance_matrix().
+   *
+   * @param[in] name1     Name of the first variable.
+   * @param[in] name2     Name of the second variable.
+   * @param[in] cov_coeff Covariance coefficient (fraction of the geometric
+   *                      mean of the two variances).
+   */
   void add_covariance_coeff(const char* name1, const char* name2,
                             double cov_coeff);
 
+  /**
+   * Return the index of a variable by name, or -1 if not found.
+   *
+   * @param[in] str Variable name to search for.
+   */
   int find_var(const std::string& str) const;
+
+  /**
+   * Update the measured value and confidence interval of an existing variable.
+   *
+   * @param[in] name    Name of the variable to update.
+   * @param[in] val     New measured value.
+   * @param[in] confint New confidence interval.
+   */
   void change_var(const char* name, double val, double confint);
+
+  /** Return the number of variables that have measurements (non-negative \c confint). */
   int get_number_measured() const;
+
+  /**
+   * Build and return the full covariance matrix from the registered variables
+   * and covariance coefficients.
+   */
   matrix get_covariance_matrix() const;
+
+  /** Return a vector of all variable values in registration order. */
   vector get_values() const;
+
+  /**
+   * Return a new \c recon_system with updated values and confidence intervals.
+   *
+   * Keeps all variable names and covariance coefficients; replaces values and
+   * confidence intervals from the supplied vectors.
+   *
+   * @param[in] values   New variable values (one per variable, in order).
+   * @param[in] confints New confidence intervals (one per variable, in order).
+   */
   recon_system updated(const vector& values, const vector& confints) const;
+
+  /** Print all variables with their values and confidence intervals to stdout. */
   void print_vars() const;
+
+  /**
+   * Print the constraint matrix \c F together with variable names to stdout.
+   *
+   * @param[in] F Constraint (Jacobian) matrix whose columns correspond to
+   *              the variables in this system.
+   */
   void print_constraints(const matrix& F) const;
 };
 
diff --git a/test/custom_reporter.h b/test/custom_reporter.h
new file mode 100644
index 0000000..e6ac338
--- /dev/null
+++ b/test/custom_reporter.h
@@ -0,0 +1,26 @@
+#pragma once
+#include <catch2/reporters/catch_reporter_event_listener.hpp>
+#include <catch2/reporters/catch_reporter_registrars.hpp>
+
+class SummaryListener : public Catch::EventListenerBase {
+public:
+  using EventListenerBase::EventListenerBase;
+
+  void testCaseEnded(Catch::TestCaseStats const& stats) override {
+    auto passed = stats.totals.assertions.passed;
+    auto failed = stats.totals.assertions.failed;
+    auto total = passed + failed;
+
+    if(failed == 0) {
+      Catch::cout() << "\033[32m  ✓ " << stats.testInfo->name
+                    << ": " << total << " assertion"
+                    << (total != 1 ? "s" : "") << " passed\033[0m\n";
+    } else {
+      Catch::cout() << "\033[31m  ✗ " << stats.testInfo->name
+                    << ": " << passed << "/" << total
+                    << " passed, " << failed << " FAILED\033[0m\n";
+    }
+  }
+};
+
+CATCH_REGISTER_LISTENER(SummaryListener)
diff --git a/test/custom_reporter.hpp b/test/custom_reporter.hpp
deleted file mode 100644
index 2444603..0000000
--- a/test/custom_reporter.hpp
+++ /dev/null
@@ -1,26 +0,0 @@
-#pragma once
-#include <catch2/reporters/catch_reporter_event_listener.hpp>
-#include <catch2/reporters/catch_reporter_registrars.hpp>
-
-class SummaryListener : public Catch::EventListenerBase {
-public:
-    using EventListenerBase::EventListenerBase;
-
-    void testCaseEnded(Catch::TestCaseStats const& stats) override {
-        auto passed = stats.totals.assertions.passed;
-        auto failed = stats.totals.assertions.failed;
-        auto total  = passed + failed;
-
-        if (failed == 0) {
-            Catch::cout() << "\033[32m  ✓ " << stats.testInfo->name
-                          << ": " << total << " assertion"
-                          << (total != 1 ? "s" : "") << " passed\033[0m\n";
-        } else {
-            Catch::cout() << "\033[31m  ✗ " << stats.testInfo->name
-                          << ": " << passed << "/" << total
-                          << " passed, " << failed << " FAILED\033[0m\n";
-        }
-    }
-};
-
-CATCH_REGISTER_LISTENER(SummaryListener)
diff --git a/test/gsl_wrapper_tests.cc b/test/gsl_wrapper_tests.cc
index 799fad7..d5b1869 100644
--- a/test/gsl_wrapper_tests.cc
+++ b/test/gsl_wrapper_tests.cc
@@ -5,218 +5,218 @@
 using namespace dvrlib;
 
 TEST_CASE("vector") {
-    vector v(3, 4.0);
-    REQUIRE(v.get(0) == 4.0);
-    REQUIRE(v.get(1) == 4.0);
-    REQUIRE(v.get(2) == 4.0);
-
-    for (int i = 0; i < 3; i++) {
-        v.set(i, 1.23 + i);
-    }
-    REQUIRE(v.get(0) == 1.23);
-    REQUIRE(v.get(1) == 2.23);
-    REQUIRE(v.get(2) == 3.23);
-
-    vector w(3, 5.0);
-    REQUIRE(w.get(0) == 5);
-
-    w += v;
-    REQUIRE(w.get(0) == 6.23);
-    REQUIRE(w.get(2) == 8.23);
-    REQUIRE(v.get(0) == 1.23);
-
-    v = w * 3;
-    REQUIRE(v.get(0) == 18.69);
-    REQUIRE(w.get(0) == 6.23);
-
-    double src[] = {3, 4, 10, 7};
-    vector x(4, src);
-    REQUIRE(x.get(0) == 3);
-    REQUIRE(x.get(3) == 7);
-
-    x -= x;
-    REQUIRE(x.get(0) == 0);
-    REQUIRE(x.get(1) == 0);
-
-    vector z(4, src);
-    REQUIRE(z.get(0) == 3);
-    REQUIRE(z.get(3) == 7);
-
-    x = -z;
-    REQUIRE(x.get(0) == -3);
-    REQUIRE(x.get(3) == -7);
-
-    REQUIRE(x[0] == -3);
-    REQUIRE(x[3] == -7);
+  vector v(3, 4.0);
+  REQUIRE(v.get(0) == 4.0);
+  REQUIRE(v.get(1) == 4.0);
+  REQUIRE(v.get(2) == 4.0);
+
+  for(int i = 0; i < 3; i++) {
+    v.set(i, 1.23 + i);
+  }
+  REQUIRE(v.get(0) == 1.23);
+  REQUIRE(v.get(1) == 2.23);
+  REQUIRE(v.get(2) == 3.23);
+
+  vector w(3, 5.0);
+  REQUIRE(w.get(0) == 5);
+
+  w += v;
+  REQUIRE(w.get(0) == 6.23);
+  REQUIRE(w.get(2) == 8.23);
+  REQUIRE(v.get(0) == 1.23);
+
+  v = w * 3;
+  REQUIRE(v.get(0) == 18.69);
+  REQUIRE(w.get(0) == 6.23);
+
+  double src[] = {3, 4, 10, 7};
+  vector x(4, src);
+  REQUIRE(x.get(0) == 3);
+  REQUIRE(x.get(3) == 7);
+
+  x -= x;
+  REQUIRE(x.get(0) == 0);
+  REQUIRE(x.get(1) == 0);
+
+  vector z(4, src);
+  REQUIRE(z.get(0) == 3);
+  REQUIRE(z.get(3) == 7);
+
+  x = -z;
+  REQUIRE(x.get(0) == -3);
+  REQUIRE(x.get(3) == -7);
+
+  REQUIRE(x[0] == -3);
+  REQUIRE(x[3] == -7);
 }
 
 TEST_CASE("vector_view") {
-    vector v(10, 4.0);
-    vector_view vv = v.subvector(3, 2);
-    vv.set(0, 2);
-    vv.set(1, 5);
-    REQUIRE(v.get(0) == 4);
-    REQUIRE(v.get(3) == 2);
-    REQUIRE(v.get(4) == 5);
-
-    vector x(vv);
-    x.set(0, 7);
-    REQUIRE(x.get(0) == 7);
-    REQUIRE(x.get(1) == 5);
-    REQUIRE(vv.get(0) == 2);
-    REQUIRE(v.get(3) == 2);
-
-    x.set(1, 8);
-    vv = x;
-    REQUIRE(v.get(3) == 7);
-    REQUIRE(v.get(4) == 8);
-
-    x = vv;
-    vv.set(0, -10);
-    REQUIRE(x.get(0) == 7);
-    REQUIRE(vv.get(0) == -10);
-
-    vector v3(3);
-    v3.set(0, 3);
-    v3.set(1, -4);
-    v3.set(2, 12);
-    REQUIRE(v3.norm1() == 19);
-    REQUIRE(v3.norm2() == 13);
+  vector v(10, 4.0);
+  vector_view vv = v.subvector(3, 2);
+  vv.set(0, 2);
+  vv.set(1, 5);
+  REQUIRE(v.get(0) == 4);
+  REQUIRE(v.get(3) == 2);
+  REQUIRE(v.get(4) == 5);
+
+  vector x(vv);
+  x.set(0, 7);
+  REQUIRE(x.get(0) == 7);
+  REQUIRE(x.get(1) == 5);
+  REQUIRE(vv.get(0) == 2);
+  REQUIRE(v.get(3) == 2);
+
+  x.set(1, 8);
+  vv = x;
+  REQUIRE(v.get(3) == 7);
+  REQUIRE(v.get(4) == 8);
+
+  x = vv;
+  vv.set(0, -10);
+  REQUIRE(x.get(0) == 7);
+  REQUIRE(vv.get(0) == -10);
+
+  vector v3(3);
+  v3.set(0, 3);
+  v3.set(1, -4);
+  v3.set(2, 12);
+  REQUIRE(v3.norm1() == 19);
+  REQUIRE(v3.norm2() == 13);
 }
 
 TEST_CASE("matrix") {
-    matrix m(5, 6, true);
-    REQUIRE(m.get(0, 0) == 1);
-    REQUIRE(m.get(0, 1) == 0);
-    REQUIRE(m.get(4, 4) == 1);
-    REQUIRE(m.get(4, 5) == 0);
-
-    m.set(3, 4, 7);
-    m.set(1, 3, 9);
-    matrix m2 = m.transpose();
-    REQUIRE(m.get(3, 4) == 7);
-    REQUIRE(m.get(1, 3) == 9);
-    REQUIRE(m2.get(4, 3) == 7);
-    REQUIRE(m2.get(3, 1) == 9);
-
-    matrix A(2, 2);
-    A.set(0, 0, 4);
-    A.set(0, 1, 6);
-    A.set(1, 0, 10);
-    A.set(1, 1, 3);
-    vector b(2);
-    b.set(0, 38);
-    b.set(1, 59);
-    vector c = A.linsolve(b);
-    REQUIRE(c.get(0) == 5);
-    REQUIRE(c.get(1) == 3);
-    REQUIRE(b.get(1) == 59);
-    REQUIRE(A.get(0, 0) == 4);
-    REQUIRE(A.get(1, 1) == 3);
-
-    matrix B(A.inverse());
-    double det = -48;
-    REQUIRE(fabs(B.get(0, 0) - +3.0 / det) < 1e-9);
-    REQUIRE(fabs(B.get(0, 1) - -6.0 / det) < 1e-9);
-    REQUIRE(fabs(B.get(1, 0) - -10.0 / det) < 1e-9);
-    REQUIRE(fabs(B.get(1, 1) - +4.0 / det) < 1e-9);
-
-    matrix A2(2, 3);
-    A2.set(0, 0, 4);
-    A2.set(0, 1, 6);
-    A2.set(0, 2, 7);
-    A2.set(1, 0, 10);
-    A2.set(1, 1, 3);
-    A2.set(1, 2, 5);
-    vector x(3);
-    x.set(0, 1);
-    x.set(1, 2);
-    x.set(2, 3);
-    vector y = A2 * x;
-    REQUIRE(y.get(0) == 37);
-    REQUIRE(y.get(1) == 31);
-
-    double src[][3] = {{3, 4, 10}, {7, 5, 11}};
-    matrix A3(2, 3, src);
-    REQUIRE(A3.get(0, 0) == 3);
-    REQUIRE(A3.get(1, 0) == 7);
-    REQUIRE(A3.get(0, 1) == 4);
-    REQUIRE(A3.get(1, 2) == 11);
-
-    matrix A4(2, 3, src);
-    A4 = A3 + A3;
-    REQUIRE(A4.get(0, 0) == 6);
-    REQUIRE(A4.get(1, 0) == 14);
-    REQUIRE(A4.get(0, 1) == 8);
-    REQUIRE(A4.get(1, 2) == 22);
-
-    A4 += A3;
-    REQUIRE(A4.get(0, 0) == 9);
-    REQUIRE(A4.get(1, 0) == 21);
-    REQUIRE(A4.get(0, 1) == 12);
-    REQUIRE(A4.get(1, 2) == 33);
-
-    A4 -= A3;
-    REQUIRE(A4.get(0, 0) == 6);
-    REQUIRE(A4.get(1, 0) == 14);
-    REQUIRE(A4.get(0, 1) == 8);
-    REQUIRE(A4.get(1, 2) == 22);
-
-    A4 = A4 - A3;
-    REQUIRE(A3.get(0, 0) == 3);
-    REQUIRE(A3.get(1, 0) == 7);
-    REQUIRE(A3.get(0, 1) == 4);
-    REQUIRE(A3.get(1, 2) == 11);
-
-    A4 = -A3;
-    REQUIRE(A3.get(0, 0) == 3);
-    REQUIRE(A3.get(1, 0) == 7);
-    REQUIRE(A3.get(0, 1) == 4);
-    REQUIRE(A3.get(1, 2) == 11);
-    REQUIRE(A4.get(0, 0) == -3);
-    REQUIRE(A4.get(1, 0) == -7);
-    REQUIRE(A4.get(0, 1) == -4);
-    REQUIRE(A4.get(1, 2) == -11);
-
-    A4 *= -2.0;
-    REQUIRE(A4.get(0, 0) == 6);
-    REQUIRE(A4.get(1, 0) == 14);
-    REQUIRE(A4.get(0, 1) == 8);
-    REQUIRE(A4.get(1, 2) == 22);
-
-    A4 = A4 * -2.0;
-    REQUIRE(A4.get(0, 0) == -12);
-    REQUIRE(A4.get(1, 0) == -28);
-    REQUIRE(A4.get(0, 1) == -16);
-    REQUIRE(A4.get(1, 2) == -44);
-
-    vector d(A4[1]);
-    REQUIRE(d.get(0) == -28);
-    REQUIRE(d.get(1) == -20);
-    REQUIRE(d.get(2) == -44);
+  matrix m(5, 6, true);
+  REQUIRE(m.get(0, 0) == 1);
+  REQUIRE(m.get(0, 1) == 0);
+  REQUIRE(m.get(4, 4) == 1);
+  REQUIRE(m.get(4, 5) == 0);
+
+  m.set(3, 4, 7);
+  m.set(1, 3, 9);
+  matrix m2 = m.transpose();
+  REQUIRE(m.get(3, 4) == 7);
+  REQUIRE(m.get(1, 3) == 9);
+  REQUIRE(m2.get(4, 3) == 7);
+  REQUIRE(m2.get(3, 1) == 9);
+
+  matrix A(2, 2);
+  A.set(0, 0, 4);
+  A.set(0, 1, 6);
+  A.set(1, 0, 10);
+  A.set(1, 1, 3);
+  vector b(2);
+  b.set(0, 38);
+  b.set(1, 59);
+  vector c = A.linsolve(b);
+  REQUIRE(c.get(0) == 5);
+  REQUIRE(c.get(1) == 3);
+  REQUIRE(b.get(1) == 59);
+  REQUIRE(A.get(0, 0) == 4);
+  REQUIRE(A.get(1, 1) == 3);
+
+  matrix B(A.inverse());
+  double det = -48;
+  REQUIRE(fabs(B.get(0, 0) - +3.0 / det) < 1e-9);
+  REQUIRE(fabs(B.get(0, 1) - -6.0 / det) < 1e-9);
+  REQUIRE(fabs(B.get(1, 0) - -10.0 / det) < 1e-9);
+  REQUIRE(fabs(B.get(1, 1) - +4.0 / det) < 1e-9);
+
+  matrix A2(2, 3);
+  A2.set(0, 0, 4);
+  A2.set(0, 1, 6);
+  A2.set(0, 2, 7);
+  A2.set(1, 0, 10);
+  A2.set(1, 1, 3);
+  A2.set(1, 2, 5);
+  vector x(3);
+  x.set(0, 1);
+  x.set(1, 2);
+  x.set(2, 3);
+  vector y = A2 * x;
+  REQUIRE(y.get(0) == 37);
+  REQUIRE(y.get(1) == 31);
+
+  double src[][3] = {{3, 4, 10}, {7, 5, 11}};
+  matrix A3(2, 3, src);
+  REQUIRE(A3.get(0, 0) == 3);
+  REQUIRE(A3.get(1, 0) == 7);
+  REQUIRE(A3.get(0, 1) == 4);
+  REQUIRE(A3.get(1, 2) == 11);
+
+  matrix A4(2, 3, src);
+  A4 = A3 + A3;
+  REQUIRE(A4.get(0, 0) == 6);
+  REQUIRE(A4.get(1, 0) == 14);
+  REQUIRE(A4.get(0, 1) == 8);
+  REQUIRE(A4.get(1, 2) == 22);
+
+  A4 += A3;
+  REQUIRE(A4.get(0, 0) == 9);
+  REQUIRE(A4.get(1, 0) == 21);
+  REQUIRE(A4.get(0, 1) == 12);
+  REQUIRE(A4.get(1, 2) == 33);
+
+  A4 -= A3;
+  REQUIRE(A4.get(0, 0) == 6);
+  REQUIRE(A4.get(1, 0) == 14);
+  REQUIRE(A4.get(0, 1) == 8);
+  REQUIRE(A4.get(1, 2) == 22);
+
+  A4 = A4 - A3;
+  REQUIRE(A3.get(0, 0) == 3);
+  REQUIRE(A3.get(1, 0) == 7);
+  REQUIRE(A3.get(0, 1) == 4);
+  REQUIRE(A3.get(1, 2) == 11);
+
+  A4 = -A3;
+  REQUIRE(A3.get(0, 0) == 3);
+  REQUIRE(A3.get(1, 0) == 7);
+  REQUIRE(A3.get(0, 1) == 4);
+  REQUIRE(A3.get(1, 2) == 11);
+  REQUIRE(A4.get(0, 0) == -3);
+  REQUIRE(A4.get(1, 0) == -7);
+  REQUIRE(A4.get(0, 1) == -4);
+  REQUIRE(A4.get(1, 2) == -11);
+
+  A4 *= -2.0;
+  REQUIRE(A4.get(0, 0) == 6);
+  REQUIRE(A4.get(1, 0) == 14);
+  REQUIRE(A4.get(0, 1) == 8);
+  REQUIRE(A4.get(1, 2) == 22);
+
+  A4 = A4 * -2.0;
+  REQUIRE(A4.get(0, 0) == -12);
+  REQUIRE(A4.get(1, 0) == -28);
+  REQUIRE(A4.get(0, 1) == -16);
+  REQUIRE(A4.get(1, 2) == -44);
+
+  vector d(A4[1]);
+  REQUIRE(d.get(0) == -28);
+  REQUIRE(d.get(1) == -20);
+  REQUIRE(d.get(2) == -44);
 }
 
 TEST_CASE("matrix_view") {
-    matrix m(10, 10, true);
-    matrix_view mv = m.submatrix(3, 3, 2, 2);
-    REQUIRE(mv.get(0, 0) == 1);
-    REQUIRE(mv.get(0, 1) == 0);
-
-    mv.set(0, 1, 2);
-    mv.set(1, 0, 5);
-    REQUIRE(mv.get(0, 1) == 2);
-    REQUIRE(mv.get(1, 0) == 5);
-    REQUIRE(m.get(3, 4) == 2);
-    REQUIRE(m.get(4, 3) == 5);
-
-    matrix m3(2, 2, true);
-    m3.set(0, 0, 5);
-    m3.set(1, 0, 6);
-    m3.set(0, 1, 7);
-    m3.set(1, 1, 8);
-    m.submatrix(2, 4, 2, 2) = m3;
-    REQUIRE(m.get(2, 4) == 5);
-    REQUIRE(m.get(3, 4) == 6);
-    REQUIRE(m.get(2, 5) == 7);
-    REQUIRE(m.get(3, 5) == 8);
+  matrix m(10, 10, true);
+  matrix_view mv = m.submatrix(3, 3, 2, 2);
+  REQUIRE(mv.get(0, 0) == 1);
+  REQUIRE(mv.get(0, 1) == 0);
+
+  mv.set(0, 1, 2);
+  mv.set(1, 0, 5);
+  REQUIRE(mv.get(0, 1) == 2);
+  REQUIRE(mv.get(1, 0) == 5);
+  REQUIRE(m.get(3, 4) == 2);
+  REQUIRE(m.get(4, 3) == 5);
+
+  matrix m3(2, 2, true);
+  m3.set(0, 0, 5);
+  m3.set(1, 0, 6);
+  m3.set(0, 1, 7);
+  m3.set(1, 1, 8);
+  m.submatrix(2, 4, 2, 2) = m3;
+  REQUIRE(m.get(2, 4) == 5);
+  REQUIRE(m.get(3, 4) == 6);
+  REQUIRE(m.get(2, 5) == 7);
+  REQUIRE(m.get(3, 5) == 8);
 }
diff --git a/test/main.cc b/test/main.cc
index 5378d85..ae67d80 100644
--- a/test/main.cc
+++ b/test/main.cc
@@ -1,10 +1,10 @@
 #include <catch2/catch_all.hpp>
 #include "gsl_wrapper.h"
-#include "custom_reporter.hpp"
+#include "custom_reporter.h"
 
 using namespace dvrlib;
 
 int main(int argc, char* argv[]) {
-    gsl_enable_exceptions();
-    return Catch::Session().run(argc, argv);
+  gsl_enable_exceptions();
+  return Catch::Session().run(argc, argv);
 }
diff --git a/test/recon_tests.cc b/test/recon_tests.cc
index 43483b6..e495449 100644
--- a/test/recon_tests.cc
+++ b/test/recon_tests.cc
@@ -6,7 +6,7 @@
 using namespace dvrlib;
 
 static bool vectors_almost_equal(const vector& a, const vector& b, double eps = 1e-7) {
-    return (a + (-1 * b)).norm2() < eps * (a.norm2() + b.norm2());
+  return (a + (-1 * b)).norm2() < eps * (a.norm2() + b.norm2());
 }
 
 namespace simple_system {
@@ -42,151 +42,158 @@ double v[] = {0.75, 0.75, -0.75, -0.75, 198.5};
 }  // namespace mischer_teiler
 
 TEST_CASE("lin_recon") {
-    SECTION("mischer_teiler") {
-        matrix F(2, 5, mischer_teiler::F);
-        matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
-        vector x(5, mischer_teiler::x);
-
-        vector r = F * x;
-        vector v(5), v_exp(5, mischer_teiler::v);
-
-        lin_recon(r, S_x, F, v);
-        REQUIRE(vectors_almost_equal(v, v_exp));
-
-        matrix S_v(S_x);
-        lin_cov_update(S_x, F, S_v);
-    }
-    SECTION("simple_system") {
-        matrix F(1, 2, simple_system::F);
-        matrix S_x(2, 2, true, simple_system::S_x_diag);
-        vector x(2, simple_system::x);
-
-        vector r = F * x;
-        vector v(2), v_exp(2, simple_system::v);
-
-        lin_recon(r, S_x, F, v);
-        REQUIRE(vectors_almost_equal(v, v_exp));
-    }
-    SECTION("simple_system2") {
-        matrix F(3, 4, simple_system2::F);
-        matrix S_x(2, 2, true, simple_system2::S_x_diag);
-        vector x(4, simple_system2::x);
-
-        vector r = F * x;
-        vector v(4), v_exp(4, simple_system2::v);
-
-        lin_recon(r, S_x, F, v);
-        REQUIRE(vectors_almost_equal(v, v_exp));
-    }
-}
-
-TEST_CASE("lin_recon_update") {
+  SECTION("mischer_teiler") {
     matrix F(2, 5, mischer_teiler::F);
-    matrix S_x_inv(4, 4, true, mischer_teiler::S_x_inv_diag);
+    matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
     vector x(5, mischer_teiler::x);
 
-    vector v(5), dv(5), v_exp(5, mischer_teiler::v);
-    v.set(1, 3);
-    v.set(2, 4);
-    v.set(4, 7);
-    vector r = F * (x + v);
+    vector r = F * x;
+    vector v(5);
+    vector v_exp(5, mischer_teiler::v);
+
+    lin_recon(r, S_x, F, v);
+    REQUIRE(vectors_almost_equal(v, v_exp));
+
+    matrix S_v(S_x);
+    lin_cov_update(S_x, F, S_v);
+  }
+  SECTION("simple_system") {
+    matrix F(1, 2, simple_system::F);
+    matrix S_x(2, 2, true, simple_system::S_x_diag);
+    vector x(2, simple_system::x);
+
+    vector r = F * x;
+    vector v(2);
+    vector v_exp(2, simple_system::v);
+
+    lin_recon(r, S_x, F, v);
+    REQUIRE(vectors_almost_equal(v, v_exp));
+  }
+  SECTION("simple_system2") {
+    matrix F(3, 4, simple_system2::F);
+    matrix S_x(2, 2, true, simple_system2::S_x_diag);
+    vector x(4, simple_system2::x);
+
+    vector r = F * x;
+    vector v(4);
+    vector v_exp(4, simple_system2::v);
+
+    lin_recon(r, S_x, F, v);
+    REQUIRE(vectors_almost_equal(v, v_exp));
+  }
+}
 
-    lin_recon_update(r, S_x_inv, F, v, dv);
-    REQUIRE(vectors_almost_equal(v + dv, v_exp));
+TEST_CASE("lin_recon_update") {
+  matrix F(2, 5, mischer_teiler::F);
+  matrix S_x_inv(4, 4, true, mischer_teiler::S_x_inv_diag);
+  vector x(5, mischer_teiler::x);
+
+  vector v(5);
+  vector v_exp(5, mischer_teiler::v);
+  vector dv(5);
+  v.set(1, 3);
+  v.set(2, 4);
+  v.set(4, 7);
+  vector r = F * (x + v);
+
+  lin_recon_update(r, S_x_inv, F, v, dv);
+  REQUIRE(vectors_almost_equal(v + dv, v_exp));
 }
 
 class EnbiproDummy {
 public:
-    virtual vector getValues() = 0;
-    virtual matrix getCovarianceMatrix() = 0;
-    virtual matrix getJacobian(const vector& x) = 0;
-    virtual vector getResidual(const vector& x) = 0;
+  virtual ~EnbiproDummy() {}
+  virtual vector getValues() = 0;
+  virtual matrix getCovarianceMatrix() = 0;
+  virtual matrix getJacobian(const vector& x) = 0;
+  virtual vector getResidual(const vector& x) = 0;
 };
 
 class LinearEnbiproDummy : public EnbiproDummy {
 public:
-    vector getValues() {
-        vector x(5, mischer_teiler::x);
-        return x;
-    }
-
-    matrix getCovarianceMatrix() {
-        matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
-        return S_x;
-    }
-
-    matrix getJacobian(const vector& x) {
-        matrix F(2, 5, mischer_teiler::F);
-        return F;
-    }
-
-    vector getResidual(const vector& x) {
-        return getJacobian(x) * x;
-    }
+  vector getValues() override {
+    vector x(5, mischer_teiler::x);
+    return x;
+  }
+
+  matrix getCovarianceMatrix() override {
+    matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
+    return S_x;
+  }
+
+  matrix getJacobian(const vector& x) override {
+    matrix F(2, 5, mischer_teiler::F);
+    return F;
+  }
+
+  vector getResidual(const vector& x) override {
+    return getJacobian(x) * x;
+  }
 };
 
 class QuadraticEnbiproDummy : public EnbiproDummy {
 public:
-    vector getValues() {
-        vector x(5, mischer_teiler::x);
-        x.set(4, 201);
-        return x;
-    }
-
-    matrix getCovarianceMatrix() {
-        matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
-        return S_x;
-    }
-
-    matrix getJacobian(const vector& x) {
-        double d1 = 2 * (x.get(0) + x.get(1) - x.get(4));
-        double d2 = 2 * (x.get(2) + x.get(3) - x.get(4));
-        double J[][5] = {{d1, d1, 0, 0, -d1}, {0, 0, d2, d2, -d2}};
-        matrix F(2, 5, J);
-        return F;
-    }
-
-    vector getResidual(const vector& x) {
-        vector r(2);
-        r.set(0, pow(x.get(0) + x.get(1) - x.get(4), 2));
-        r.set(1, pow(-x.get(2) - x.get(3) + x.get(4), 2));
-        return r;
-    }
+  vector getValues() override {
+    vector x(5, mischer_teiler::x);
+    x.set(4, 201);
+    return x;
+  }
+
+  matrix getCovarianceMatrix() override {
+    matrix S_x(4, 4, true, mischer_teiler::S_x_diag);
+    return S_x;
+  }
+
+  matrix getJacobian(const vector& x) override {
+    double d1 = 2 * (x.get(0) + x.get(1) - x.get(4));
+    double d2 = 2 * (x.get(2) + x.get(3) - x.get(4));
+    double J[][5] = {{d1, d1, 0, 0, -d1}, {0, 0, d2, d2, -d2}};
+    matrix F(2, 5, J);
+    return F;
+  }
+
+  vector getResidual(const vector& x) override {
+    vector r(2);
+    r.set(0, pow(x.get(0) + x.get(1) - x.get(4), 2));
+    r.set(1, pow(-x.get(2) - x.get(3) + x.get(4), 2));
+    return r;
+  }
 };
 
 class EnbiJacobianFunc : public func<vector, matrix> {
-    EnbiproDummy* enbi;
+  EnbiproDummy* enbi;
 
 public:
-    EnbiJacobianFunc(EnbiproDummy* _enbi) : enbi(_enbi) {}
-    virtual matrix operator()(const vector& arg) const;
+  EnbiJacobianFunc(EnbiproDummy* _enbi) : enbi(_enbi) {}
+  matrix operator()(const vector& arg) const override;
 };
 matrix EnbiJacobianFunc::operator()(const vector& arg) const {
-    return enbi->getJacobian(arg);
+  return enbi->getJacobian(arg);
 }
 
 class EnbiResidualFunc : public func<vector, vector> {
-    EnbiproDummy* enbi;
+  EnbiproDummy* enbi;
 
 public:
-    EnbiResidualFunc(EnbiproDummy* _enbi) : enbi(_enbi) {}
-    virtual vector operator()(const vector& arg) const;
+  EnbiResidualFunc(EnbiproDummy* _enbi) : enbi(_enbi) {}
+  vector operator()(const vector& arg) const override;
 };
 vector EnbiResidualFunc::operator()(const vector& arg) const {
-    return enbi->getResidual(arg);
+  return enbi->getResidual(arg);
 }
 
 TEST_CASE("recon nonlinear") {
-    QuadraticEnbiproDummy enbi_dummy;
-    vector x = enbi_dummy.getValues();
-    matrix S_x = enbi_dummy.getCovarianceMatrix();
-    EnbiResidualFunc f(&enbi_dummy);
-    EnbiJacobianFunc J(&enbi_dummy);
-    vector v(x.size()), v_exp(5, mischer_teiler::v);
-    v_exp.set(4, -2.5);  // note that x(4) is different here, thus v(4) also
-    matrix S_v(v.size(), v.size());
-
-    recon(x, S_x, f, J, v, S_v);
-    REQUIRE(f(x + v).norm2() < 1e-6);
-    REQUIRE(vectors_almost_equal(v, v_exp, 1e-3));
+  QuadraticEnbiproDummy enbi_dummy;
+  vector x = enbi_dummy.getValues();
+  matrix S_x = enbi_dummy.getCovarianceMatrix();
+  EnbiResidualFunc f(&enbi_dummy);
+  EnbiJacobianFunc J(&enbi_dummy);
+  vector v(x.size());
+  vector v_exp(5, mischer_teiler::v);
+  v_exp.set(4, -2.5);  // note that x(4) is different here, thus v(4) also
+  matrix S_v(v.size(), v.size());
+
+  recon(x, S_x, f, J, v, S_v);
+  REQUIRE(f(x + v).norm2() < 1e-6);
+  REQUIRE(vectors_almost_equal(v, v_exp, 1e-3));
 }
diff --git a/test/vdi2048_test.cc b/test/vdi2048_test.cc
index e732de6..50d0600 100644
--- a/test/vdi2048_test.cc
+++ b/test/vdi2048_test.cc
@@ -6,102 +6,99 @@
 using namespace dvrlib;
 
 static void setup_system(recon_system& sys) {
-    sys.add_var("m_FDKeI",  46.241, 0.800);
-    sys.add_var("m_FDKeII", 45.668, 0.790);
-    sys.add_var("m_SpI",    44.575, 0.535);
-    sys.add_var("m_SpII",   44.319, 0.532);
-    sys.add_var("m_V",       0.525, 0.105);
-    sys.add_var("m_HK",     69.978, 0.854);
-    sys.add_var("m_A7",     10.364, 0.168);
-    sys.add_var("m_A6",      3.744, 0.058);
-    sys.add_var("m_A5",      4.391, 0.058);
-    sys.add_var("m_HDNK",   18.498, 0.205);
-    sys.add_var("m_D",       2.092, 0.272);
-    sys.add_var("m_FD1",     0,    -1);
-    sys.add_var("m_FD2",     0,    -1);
-    sys.add_var("m_FD3",     0,    -1);
-    sys.add_var("m_HDAnz",   0,    -1);
-    sys.add_covariance_coeff("m_FDKeI",  "m_FDKeII", 0.2);
-    sys.add_covariance_coeff("m_SpI",    "m_SpII",   0.4);
+  sys.add_var("m_FDKeI",  46.241, 0.800);
+  sys.add_var("m_FDKeII", 45.668, 0.790);
+  sys.add_var("m_SpI",    44.575, 0.535);
+  sys.add_var("m_SpII",   44.319, 0.532);
+  sys.add_var("m_V",       0.525, 0.105);
+  sys.add_var("m_HK",     69.978, 0.854);
+  sys.add_var("m_A7",     10.364, 0.168);
+  sys.add_var("m_A6",      3.744, 0.058);
+  sys.add_var("m_A5",      4.391, 0.058);
+  sys.add_var("m_HDNK",   18.498, 0.205);
+  sys.add_var("m_D",       2.092, 0.272);
+  sys.add_var("m_FD1",     0,  -1);
+  sys.add_var("m_FD2",     0,  -1);
+  sys.add_var("m_FD3",     0,  -1);
+  sys.add_var("m_HDAnz",   0,  -1);
+  sys.add_covariance_coeff("m_FDKeI",  "m_FDKeII", 0.2);
+  sys.add_covariance_coeff("m_SpI",  "m_SpII",   0.4);
 }
 
-static matrix make_F() {
-    double Fc[][15] = {
-        {1, 1, 0, 0, -0.2, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0},
-        {0, 0, 1, 1, -0.6, 0, 0, 0, 0, 0, 0,  0,-1, 0, 0},
-        {0, 0, 0, 0,  0.4, 1, 1, 1, 1, 0, 0,  0, 0,-1, 0},
-        {0, 0, 0, 0,  0,   0, 1, 1, 1, 0, 0,  0, 0, 0,-1},
-        {0, 0, 0, 0,  0,   0, 0, 0, 0, 0, 0,  1,-1, 0, 0},
-        {0, 0, 0, 0,  0,   0, 0, 0, 0, 0, 0,  0, 1,-1, 0},
-        {0, 0, 0, 0,  0,   0, 0, 0, 0,-1, 0,  0, 0, 0, 1},
-    };
-    return matrix(7, 15, Fc);
+static matrix make_constraints_matrix() {
+  double Fc[][15] = {
+    {1, 1, 0, 0, -0.2, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0},
+    {0, 0, 1, 1, -0.6, 0, 0, 0, 0, 0, 0,  0,-1, 0, 0},
+    {0, 0, 0, 0,  0.4, 1, 1, 1, 1, 0, 0,  0, 0,-1, 0},
+    {0, 0, 0, 0,  0,   0, 1, 1, 1, 0, 0,  0, 0, 0,-1},
+    {0, 0, 0, 0,  0,   0, 0, 0, 0, 0, 0,  1,-1, 0, 0},
+    {0, 0, 0, 0,  0,   0, 0, 0, 0, 0, 0,  0, 1,-1, 0},
+    {0, 0, 0, 0,  0,   0, 0, 0, 0,-1, 0,  0, 0, 0, 1},
+  };
+  return {7, 15, Fc};
 }
 
 TEST_CASE("VDI2048 regression") {
-    recon_system sys;
-    setup_system(sys);
-    matrix S_x = sys.get_covariance_matrix();
-    matrix F   = make_F();
-    vector x   = sys.get_values();
+  recon_system sys;
+  setup_system(sys);
+  matrix S_x = sys.get_covariance_matrix();
+  matrix F = make_constraints_matrix();
+  vector x = sys.get_values();
 
-    vector v(x.size());
-    lin_recon(F * x, S_x, F, v);
+  vector v(x.size());
+  lin_recon(F * x, S_x, F, v);
 
-    matrix S_v(S_x.size1(), S_x.size2());
-    lin_cov_update(S_x, F, S_v);
+  matrix S_v(S_x.size1(), S_x.size2());
+  lin_cov_update(S_x, F, S_v);
 
-    matrix S_x_new = S_x - S_v;
-    vector x_new   = x + v;
-    vector confint_new(S_x.size1());
-    extract_confidence(S_x_new, confint_new);
+  matrix S_x_new = S_x - S_v;
+  vector x_new = x + v;
+  vector confint_new(S_x.size1());
+  extract_confidence(S_x_new, confint_new);
 
-    const double tol = 1e-6;
+  constexpr auto tol = 1e-6;
 
-    SECTION("corrections v") {
-        double expected[] = {
-            -1.27651079357085, -1.25002785826582,
-             0.350904047191697, 0.348096449940544,
-             0.00134787242233519, 0.575060829315,
-             0.0117829833264503, 0.00140440603423252,
-             0.00140440603423252, 0.0155917953949165,
-             0, 89.2771917736789, 89.2771917736789,
-             89.2771917736788, 18.5135917953949
-        };
-        for (int i = 0; i < v.size(); i++)
-            CHECK(v.get(i) == Catch::Approx(expected[i]).epsilon(tol));
-    }
+  SECTION("corrections v") {
+    double expected[] = {
+      -1.27651079357085, -1.25002785826582,
+      0.350904047191697, 0.348096449940544,
+      0.00134787242233519, 0.575060829315,
+      0.0117829833264503, 0.00140440603423252,
+      0.00140440603423252, 0.0155917953949165,
+      0, 89.2771917736789, 89.2771917736789,
+      89.2771917736788, 18.5135917953949};
+    for(int i = 0; i < v.size(); i++)
+      CHECK(v.get(i) == Catch::Approx(expected[i]).epsilon(tol));
+  }
 
-    SECTION("constraints satisfied after correction") {
-        vector residual = F * (x + v);
-        for (int i = 0; i < residual.size(); i++)
-            CHECK(residual.get(i) == Catch::Approx(0.0).margin(1e-10));
-    }
+  SECTION("constraints satisfied after correction") {
+    vector residual = F * (x + v);
+    for(int i = 0; i < residual.size(); i++)
+      CHECK(residual.get(i) == Catch::Approx(0.0).margin(1e-10));
+  }
 
-    SECTION("updated values x_new") {
-        double expected[] = {
-            44.9644892064291, 44.4179721417342,
-            44.9259040471917, 44.6670964499405,
-             0.526347872422335, 70.553060829315,
-            10.3757829833265,  3.74540440603423,
-             4.39240440603423, 18.5135917953949,
-             2.092, 89.2771917736789, 89.2771917736789,
-            89.2771917736788, 18.5135917953949
-        };
-        for (int i = 0; i < x_new.size(); i++)
-            CHECK(x_new.get(i) == Catch::Approx(expected[i]).epsilon(tol));
-    }
+  SECTION("updated values x_new") {
+    double expected[] = {
+      44.9644892064291, 44.4179721417342,
+      44.9259040471917, 44.6670964499405,
+      0.526347872422335, 70.553060829315,
+      10.3757829833265, 3.74540440603423,
+      4.39240440603423, 18.5135917953949,
+      2.092, 89.2771917736789, 89.2771917736789,
+      89.2771917736788, 18.5135917953949};
+    for(int i = 0; i < x_new.size(); i++)
+      CHECK(x_new.get(i) == Catch::Approx(expected[i]).epsilon(tol));
+  }
 
-    SECTION("updated confidence intervals") {
-        double expected[] = {
-            0.575626570589129, 0.572817587480311,
-            0.404463177421466, 0.402919370199657,
-            0.104623410451298, 0.559872848891984,
-            0.133003389109488, 0.056695251568643,
-            0.056695251568643, 0.137102511396334,
-            0.272
-        };
-        for (int i = 0; i < confint_new.size(); i++)
-            CHECK(confint_new.get(i) == Catch::Approx(expected[i]).epsilon(tol));
-    }
+  SECTION("updated confidence intervals") {
+    double expected[] = {
+      0.575626570589129, 0.572817587480311,
+      0.404463177421466, 0.402919370199657,
+      0.104623410451298, 0.559872848891984,
+      0.133003389109488, 0.056695251568643,
+      0.056695251568643, 0.137102511396334,
+      0.272};
+    for(int i = 0; i < confint_new.size(); i++)
+      CHECK(confint_new.get(i) == Catch::Approx(expected[i]).epsilon(tol));
+  }
 }