diff --git a/build/Makefile.am b/build/Makefile.am
index 77eed4835f4de57ed09021fcc6f581b4792cd9b2..6b54b8e6315e8b55bc1ad3d882bdca28f205eb6e 100644
--- a/build/Makefile.am
+++ b/build/Makefile.am
@@ -1,6 +1,6 @@
 LDADD=libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}
 lib_LTLIBRARIES=libnptm/libnptm.la
-libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
+libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
 if BUILDFORTRAN
 PROGS=cluster/edfb_clu cluster/clu cluster/np_cluster sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS
 bin_PROGRAMS=$(PROGS)
diff --git a/build/Makefile.in b/build/Makefile.in
index ac105e2334fff5ff736df9e6914a58a2736f72a1..dfc4f1b8f9bb056082ad149aa35ba966718afbcd 100644
--- a/build/Makefile.in
+++ b/build/Makefile.in
@@ -156,8 +156,9 @@ am_libnptm_libnptm_la_OBJECTS = ../src/libnptm/algebraic.lo \
 	../src/libnptm/Configuration.lo ../src/libnptm/file_io.lo \
 	../src/libnptm/lapack_calls.lo ../src/libnptm/logging.lo \
 	../src/libnptm/magma_calls.lo ../src/libnptm/Parsers.lo \
-	../src/libnptm/sph_subs.lo ../src/libnptm/tfrfme.lo \
-	../src/libnptm/TransitionMatrix.lo ../src/libnptm/tra_subs.lo
+	../src/libnptm/sph_subs.lo ../src/libnptm/utils.lo \
+	../src/libnptm/tfrfme.lo ../src/libnptm/TransitionMatrix.lo \
+	../src/libnptm/tra_subs.lo
 libnptm_libnptm_la_OBJECTS = $(am_libnptm_libnptm_la_OBJECTS)
 AM_V_lt = $(am__v_lt_@AM_V@)
 am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@)
@@ -304,6 +305,7 @@ am__depfiles_remade = ../src/cluster/$(DEPDIR)/cluster.Po \
 	../src/libnptm/$(DEPDIR)/sph_subs.Plo \
 	../src/libnptm/$(DEPDIR)/tfrfme.Plo \
 	../src/libnptm/$(DEPDIR)/tra_subs.Plo \
+	../src/libnptm/$(DEPDIR)/utils.Plo \
 	../src/sphere/$(DEPDIR)/np_sphere.Po \
 	../src/sphere/$(DEPDIR)/sphere.Po \
 	../src/testing/$(DEPDIR)/test_TEDF.Po \
@@ -454,6 +456,7 @@ CXXCPP = @CXXCPP@
 CXXDEPMODE = @CXXDEPMODE@
 CXXFLAGS = @CXXFLAGS@
 CYGPATH_W = @CYGPATH_W@
+DEBUGFLAGS = @DEBUGFLAGS@
 DEFS = @DEFS@
 DEPDIR = @DEPDIR@
 DLLTOOL = @DLLTOOL@
@@ -577,7 +580,7 @@ top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 LDADD = libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}
 lib_LTLIBRARIES = libnptm/libnptm.la
-libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
+libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp
 @BUILDFORTRAN_FALSE@PROGS = cluster/np_cluster sphere/np_sphere trapping/np_trapping testing/test_TEDF testing/test_TTMS
 @BUILDFORTRAN_TRUE@PROGS = cluster/edfb_clu cluster/clu cluster/np_cluster sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_TEDF testing/test_TTMS
 @BUILDFORTRAN_TRUE@EDFBCLUSOURCES = ../src/cluster/edfb_clu.f
@@ -749,6 +752,8 @@ clean-libLTLIBRARIES:
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/sph_subs.lo: ../src/libnptm/$(am__dirstamp) \
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
+../src/libnptm/utils.lo: ../src/libnptm/$(am__dirstamp) \
+	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/tfrfme.lo: ../src/libnptm/$(am__dirstamp) \
 	../src/libnptm/$(DEPDIR)/$(am__dirstamp)
 ../src/libnptm/TransitionMatrix.lo: ../src/libnptm/$(am__dirstamp) \
@@ -900,6 +905,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/sph_subs.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/tfrfme.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/tra_subs.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/utils.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/sphere/$(DEPDIR)/np_sphere.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/sphere/$(DEPDIR)/sphere.Po@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@../src/testing/$(DEPDIR)/test_TEDF.Po@am__quote@ # am--include-marker
@@ -1268,6 +1274,7 @@ distclean: distclean-am
 	-rm -f ../src/libnptm/$(DEPDIR)/sph_subs.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tfrfme.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tra_subs.Plo
+	-rm -f ../src/libnptm/$(DEPDIR)/utils.Plo
 	-rm -f ../src/sphere/$(DEPDIR)/np_sphere.Po
 	-rm -f ../src/sphere/$(DEPDIR)/sphere.Po
 	-rm -f ../src/testing/$(DEPDIR)/test_TEDF.Po
@@ -1337,6 +1344,7 @@ maintainer-clean: maintainer-clean-am
 	-rm -f ../src/libnptm/$(DEPDIR)/sph_subs.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tfrfme.Plo
 	-rm -f ../src/libnptm/$(DEPDIR)/tra_subs.Plo
+	-rm -f ../src/libnptm/$(DEPDIR)/utils.Plo
 	-rm -f ../src/sphere/$(DEPDIR)/np_sphere.Po
 	-rm -f ../src/sphere/$(DEPDIR)/sphere.Po
 	-rm -f ../src/testing/$(DEPDIR)/test_TEDF.Po
diff --git a/build/configure b/build/configure
index 5fa0692ffb1fd46bed6a5ab4ad5bd3d331f365a9..72b47b22d9fe3727ffb8a279f2c62efa15d53fae 100755
--- a/build/configure
+++ b/build/configure
@@ -669,6 +669,7 @@ OFFLOADFLAGS
 BUILDFORTRAN_FALSE
 BUILDFORTRAN_TRUE
 ENABLE_ILP64
+DEBUGFLAGS
 HDF5_LDFLAGS
 CXXCPP
 LT_SYS_LIBRARY_PATH
@@ -813,6 +814,7 @@ with_aix_soname
 with_gnu_ld
 with_sysroot
 enable_libtool_lock
+enable_debug
 enable_ilp64
 enable_fortran
 enable_offload
@@ -1479,6 +1481,7 @@ Optional Features:
   --enable-fast-install[=PKGS]
                           optimize for fast installation [default=yes]
   --disable-libtool-lock  avoid locking (might break parallel builds)
+  --enable-debug=FEATURE  enable debug of FEATURE [default=no]
   --enable-ilp64          enable 64-bit indexing [default=yes]
   --enable-fortran        enable legacy FORTRAN compilation [default=auto]
   --enable-offload        enable target offloading (requires g++ version >=
@@ -24936,6 +24939,25 @@ esac
 fi
 
 # Configure the optional features
+# Check whether --enable-debug was given.
+if test ${enable_debug+y}
+then :
+  enableval=$enable_debug;
+    if test "x$enableval" = "xAM"; then
+      DEBUGFLAGS="-DDEBUG_AM"
+
+    fi
+
+else case e in #(
+  e)
+    DEBUGFLAGS=""
+
+
+ ;;
+esac
+fi
+
+
 # Check whether --enable-ilp64 was given.
 if test ${enable_ilp64+y}
 then :
@@ -25647,7 +25669,7 @@ else case e in #(
  ;;
 esac
 fi
-CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
+CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $DEBUGFLAGS $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
 SUBDIRS="cluster libnptm sphere testing trapping"
 
 # Generate the output
diff --git a/build/configure.ac b/build/configure.ac
index 7ff1a35c695fef8c5c07baefcd69ee239d4e48eb..0cc9b84493d563d6f9466f27a2bf1d1b01612cbb 100644
--- a/build/configure.ac
+++ b/build/configure.ac
@@ -405,6 +405,19 @@ AS_IF(
 )
 
 # Configure the optional features
+AC_ARG_ENABLE(
+  [debug],
+  [AS_HELP_STRING([--enable-debug=FEATURE], [enable debug of FEATURE [default=no]])],
+  [
+    if test "x$enableval" = "xAM"; then
+      AC_SUBST([DEBUGFLAGS], ["-DDEBUG_AM"])
+    fi
+  ],
+  [
+    AC_SUBST([DEBUGFLAGS], [""])
+  ]
+)
+
 AC_ARG_ENABLE(
   [ilp64],
   [AS_HELP_STRING([--enable-ilp64], [enable 64-bit indexing [default=yes]])],
@@ -641,7 +654,7 @@ AS_IF(
   [AC_SUBST([OMPFLAGS], [""])],
   [AC_SUBST([OMPFLAGS], [$OMPFLAGS])]
 )
-CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
+CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $DEBUGFLAGS $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
 SUBDIRS="cluster libnptm sphere testing trapping"
 
 # Generate the output
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 6a2f841fba6b496db255caa32bbc7f47a3a7457b..e7da6b93ad2f1bd4fb51d352bfc6195beecf3d28 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -83,6 +83,10 @@
 #include "../include/file_io.h"
 #endif
 
+#ifndef INCLUDE_UTILS_H_
+#include "../include/utils.h"
+#endif
+
 using namespace std;
 
 // I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
@@ -370,7 +374,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
       int myjxi488startoffset = 0;
       int myMPIstride = ompnumthreads;
       int myMPIblock = ompnumthreads;
-      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access the all later
+      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later
       VirtualAsciiFile **p_outarray = NULL;
       VirtualBinaryFile **vtppoanarray = NULL;
 
@@ -762,8 +766,31 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
   nvtxRangePush("Calculate inverted matrix");
+#endif
+#ifdef DEBUG_AM
+  /* now, before cms, output am to p_outam0 */
+  VirtualAsciiFile *outam0 = new VirtualAsciiFile();
+  string outam0_name = output_path + "/c_AM0_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix before CMS\n");
+  outam0->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam0->append_line(virtual_line);
+  write_dcomplex_matrix(outam0, cid->am, ndit, ndit);
+  outam0->write_to_disk(outam0_name);
+  delete outam0;
 #endif
   cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam1 = new VirtualAsciiFile();
+  string outam1_name = output_path + "/c_AM1_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after CMS before LUCIN\n");
+  outam1->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam1->append_line(virtual_line);
+  write_dcomplex_matrix(outam1, cid->am, ndit, ndit, " %5d %5d (%17.8lE,%17.8lE)\n", 1);
+  outam1->write_to_disk(outam1_name);
+  delete outam1;
+#endif
 #ifdef USE_NVTX
   nvtxRangePop();
 #endif
@@ -776,6 +803,17 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   nvtxRangePush("Invert the matrix");
 #endif
   invert_matrix(cid->am, ndit, jer, mxndm, cid->proc_device);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam2 = new VirtualAsciiFile();
+  string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after LUCIN before ZTM\n");
+  outam2->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam2->append_line(virtual_line);
+  write_dcomplex_matrix(outam2, cid->am, ndit, ndit);
+  outam2->write_to_disk(outam2_name);
+  delete outam2;
+#endif
 #ifdef USE_NVTX
   nvtxRangePop();
 #endif
@@ -794,6 +832,17 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   nvtxRangePush("Average calculation");
 #endif
   ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam3 = new VirtualAsciiFile();
+  string outam3_name = output_path + "/c_AM3_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after ZTM\n");
+  outam3->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam3->append_line(virtual_line);
+  write_dcomplex_matrix(outam3, cid->am, ndit, ndit);
+  outam3->write_to_disk(outam3_name);
+  delete outam3;
+#endif
   if (idfc >= 0) {
     if (jxi488 == jwtm) {
       int nlemt = 2 * cid->c4->nlem;
diff --git a/src/include/utils.h b/src/include/utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..82779beb33fa14026b8a20c3cbf9f47a01a0c65c
--- /dev/null
+++ b/src/include/utils.h
@@ -0,0 +1,40 @@
+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.h
+ *
+ * \brief Definition of auxiliary code utilities.
+ */
+
+#ifndef INCLUDE_UTILS_H_
+#define INCLUDE_UTILS_H_
+
+/*! \brief Write a double complex matrix to a text file.
+ *
+ * \param af: `VirtualAsciiFile *` Pointer to an existing VirtualAsciiFile.
+ * \param mat: `dcomplex **` Pointer to the matrix.
+ * \param rows: `int` Number of rows in the matrix.
+ * \param columns: `int` Number of columns in the matrix.
+ * \param format: `const string&` Format of the line (default is \" %5d %5d (%17.8lE,%17.8lE)\n\")
+ * \param first_index: `int` Index of the first element (default is 1, i.e. base 1 FORTRAN array notation)
+ */
+int write_dcomplex_matrix(
+			  VirtualAsciiFile *af, dcomplex **mat, int rows,
+			  int columns, const std::string& format=" %5d %5d (%17.8lE,%17.8lE)\n",
+			  int first_index=1
+);
+
+#endif
diff --git a/src/libnptm/clu_subs.cpp b/src/libnptm/clu_subs.cpp
index 42d136b28844f9ddb2a13f1909e355d1e742918e..78bfc4af3adbbc82f73387b8125347f34e6eaacd 100644
--- a/src/libnptm/clu_subs.cpp
+++ b/src/libnptm/clu_subs.cpp
@@ -1473,21 +1473,25 @@ void polar(
   bool onx = (y == 0.0);
   bool ony = (x == 0.0);
   bool onz = (onx && ony);
+  double rhos = 0.0;
   double rho = 0.0;
   if (!onz) {
     if (!onx) {
       if (!ony) {
-	rho = sqrt(x * x + y * y);
+	rhos = x * x + y * y;
+	rho = sqrt(rhos);
 	cph = x / rho;
 	sph = y / rho;
 	// goes to 25
       } else { // label 20
+	rhos = y * y;
 	rho = (y > 0.0) ? y : -y;
 	cph = 0.0;
 	sph = (y > 0.0) ? 1.0 : -1.0;
 	// goes to 25
       }
     } else { // label 15
+      rhos = x * x;
       rho = (x > 0.0) ? x : -x;
       cph = (x > 0.0) ? 1.0 : -1.0;
       sph = 0.0;
@@ -1513,7 +1517,7 @@ void polar(
     }
   } else { // label 35
     if (!onz) {
-      r = sqrt(rho * rho + z * z);
+      r = sqrt(rhos + z * z);
       cth = z / r;
       sth = rho / r;
       // returns
diff --git a/src/libnptm/sph_subs.cpp b/src/libnptm/sph_subs.cpp
index 4a63948008c71a3362bc6b25d44818a25e7efa0f..046c280d466b5b521607b1d89713a3670f76b82a 100644
--- a/src/libnptm/sph_subs.cpp
+++ b/src/libnptm/sph_subs.cpp
@@ -123,7 +123,7 @@ void cbf(int n, dcomplex z, int &nm, dcomplex *csj) {
   if (m < n) nm = m;
   else m = msta2(a0, n, 15);
   dcomplex cf0 = 0.0;
-  dcomplex cf1 = 1.0e-100;
+  dcomplex cf1 = 1.0e+00 - 100.0;
   dcomplex cf, cs;
   for (int k = m; k >= 0; k--) {
     cf = (2.0 * k + 3.0) * cf1 / z - cf0;
@@ -668,7 +668,7 @@ void rbf(int n, double x, int &nm, double sj[]) {
   if (m < n) nm = m;
   else m = msta2(a0, n, 15);
   double f0 = 0.0;
-  double f1 = 1.0e-100;
+  double f1 = 1.0e+00 - 100.0;
   double f;
   for (int k = m; k >= 0; k--) {
     f = (2.0 * k +3.0) * f1 / x - f0;
diff --git a/src/libnptm/utils.cpp b/src/libnptm/utils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f989f07565b33e18d93ed6390675468a63b41fb3
--- /dev/null
+++ b/src/libnptm/utils.cpp
@@ -0,0 +1,82 @@
+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.cpp
+ *
+ * \brief Implementation of auxiliary code utilities.
+ */
+
+/* Copyright (C) 2024   INAF - Osservatorio Astronomico di Cagliari
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+   
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+   
+   A copy of the GNU General Public License is distributed along with
+   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.
+ */
+
+/*! \file utils.h
+ *
+ * \brief Definition of auxiliary code utilities.
+ */
+
+#include <hdf5.h>
+#ifndef INCLUDE_TYPES_H_
+#include "../include/types.h"
+#endif
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
+#ifndef INCLUDE_UTILS_H_
+#include "../include/utils.h"
+#endif
+
+using namespace std;
+
+int write_dcomplex_matrix(
+			  VirtualAsciiFile *af, dcomplex **mat, int rows, int columns,
+			  const std::string& format, int first_index
+) {
+  int result = 0;
+  char virtual_line[256];
+  for (int i=0; i < rows; i++) {
+    for (int j = 0; j < columns; j++) {
+      sprintf(
+	      virtual_line, format.c_str(), i + first_index, j + first_index,
+	      real(mat[i][j]), imag(mat[i][j])
+      );
+      af->append_line(virtual_line);
+    }
+  }
+  return result;
+}