diff --git a/src/include/clu_subs.h b/src/include/clu_subs.h
index 9bc2f72a9c50808d78d9f035fac89cf22446f217..49b89e20875919f916cea21a269fa7b210da3a9f 100644
--- a/src/include/clu_subs.h
+++ b/src/include/clu_subs.h
@@ -15,7 +15,11 @@
#ifndef INCLUDE_CLU_SUBS_H_
#define INCLUDE_CLU_SUBS_H_
-/*! \brief C++ porting of APC
+/*! \brief Compute the asymmetry-corrected scattering cross-section of a cluster.
+ *
+ * This function computes the product between the geometrical asymmetry parameter and
+ * the scattering cross-section, like `aps()`, but for a cluster of spheres. See Eq. (3.16)
+ * in Borghese, Denti & Saija (2007).
*
* \param zpv: `double ****`
* \param le: `int`
@@ -30,7 +34,11 @@ void apc(
double sqk, double **gapr, std::complex<double> **gapp
);
-/*! \brief C++ porting of APCRA
+/*! \brief Compute the asymmetry-corrected scattering cross-section under random average
+ * conditions.
+ *
+ * This function computes the product between the geometrical asymmetry parameter and
+ * the scattering cross-section of a cluster using the random average directions.
*
* \param zpv: `double ****`
* \param le: `int`
@@ -45,7 +53,9 @@ void apcra(
double **gaprm, std::complex<double> **gappm
);
-/*! \brief C++ porting of CDTP
+/*! \brief Complex inner product.
+ *
+ * This function performs the complex inner product. It is used by `lucin()`.
*
* \param z: `complex<double>`
* \param am: Matrix of complex.
@@ -59,7 +69,7 @@ std::complex<double> cdtp(
int k, int nj
);
-/*! \brief C++ porting of CGEV
+/*! \brief C++ porting of CGEV. QUESTION: description?
*
* \param ipamo: `int`
* \param mu: `int`
@@ -69,7 +79,10 @@ std::complex<double> cdtp(
*/
double cgev(int ipamo, int mu, int l, int m);
-/*! \brief C++ porting of CMS
+/*! \brief Build the multi-centered M-matrix of the cluster.
+ *
+ * This function constructs the multi-centered M-matrix of the cluster, according
+ * to Eq. (5.28) of Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param c1: `C1 *`
@@ -79,7 +92,11 @@ double cgev(int ipamo, int mu, int l, int m);
*/
void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
-/*! \brief C++ porting of CRSM1
+/*! \brief Compute orientation-averaged scattered field intensity.
+ *
+ * This function computes the intensity of the scattered field for the cluster,
+ * averaged on the orientations. It is invoked for IAVM=1 (geometry referred to
+ * the meridional plane). QUESTION: correct?
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -90,7 +107,10 @@ void cms(std::complex<double> **am, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
*/
void crsm1(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4, C6 *c6);
-/*! \brief C++ porting of GHIT
+/*! \brief Compute the transfer vector from N2 to N1.
+ *
+ * This function computes the transfer vector going from N2 to N1, using either
+ * Hankel, Bessel or Bessel from origin functions.
*
* \param ihi: `int`
* \param ipamo: `int`
@@ -109,7 +129,10 @@ std::complex<double> ghit(
C1_AddOns *c1ao, C4 *c4, C6 *c6
);
-/*! \brief C++ porting of HJV
+/*! \brief Compute Hankel funtion and Bessel functions.
+ *
+ * This function constructs the Hankel function and the Bessel functions vectors. See
+ * page 331 in Borghese, Denti & Saija (2007).
*
* \param exri: `double` External medium refractive index.
* \param vk: `double` Wave number.
@@ -125,7 +148,10 @@ void hjv(
C1 *c1, C1_AddOns *c1ao, C4 *c4
);
-/*! \brief C++ porting of LUCIN
+/*! \brief Invert the multi-centered M-matrix.
+ *
+ * This function performs the inversion of the multi-centered M-matrix through
+ * LU decomposition. See Eq. (5.29) in Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param nddmst: `const int`
@@ -134,7 +160,11 @@ void hjv(
*/
void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
-/*! \brief C++ porting of MEXTC
+/*! \brief Compute the average extinction cross-section.
+ *
+ * This funbction computes the average extinction cross-section starting
+ * from the definition of the scattering amplitude. See Sec. 3.2.1 of Borghese,
+ * Denti & Saija (2007).
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -144,7 +174,10 @@ void lucin(std::complex<double> **am, const int nddmst, int n, int &ier);
*/
void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr, double **cext);
-/*! \brief C++ porting of PCROS
+/*! \brief Compute cross-sections and forward scattering amplitude for the cluster.
+ *
+ * This function computes the scattering, absorption and extinction cross-sections
+ * of the cluster, together with the Forward Scattering Amplitude.
*
* This function is intended to evaluate the particle cross-section. QUESTIUON: correct?
* \param vk: `double` Wave number.
@@ -155,7 +188,10 @@ void mextc(double vk, double exri, std::complex<double> **fsac, double **cextlr,
*/
void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4);
-/*! \brief C++ porting of PCRSM0
+/*! \brief Compute orientation-averaged cross-sections and forward scattering amplitude.
+ *
+ * This function computes the orientation-averaged scattering, absorption and extinction
+ * cross-sections of the cluster, together with the averaged Forward Scattering Amplitude.
*
* \param vk: `double` Wave number.
* \param exri: `double` External medium refractive index.
@@ -166,7 +202,10 @@ void pcros(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C4 *c4);
*/
void pcrsm0(double vk, double exri, int inpol, C1 *c1, C1_AddOns *c1ao, C4 *c4);
-/*! \brief C++ porting of POLAR
+/*! \brief Transform Cartesian quantities to spherical coordinates ones.
+ *
+ * This function performs a conversion from the Cartesian coordinates system to the
+ * spherical one. It is used by `sphar()`.
*
* \param x: `double` X-axis Cartesian coordinate.
* \param y: `double` Y-axis Cartesian coordinate.
@@ -182,7 +221,10 @@ void polar(
double &cph, double &sph
);
-/*! \brief C++ porting of R3J000
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients towards J=0.
+ *
+ * This function calculates the 3j(J,J2,J3;0,0,0) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j2: `int`
* \param j3: `int`
@@ -190,7 +232,10 @@ void polar(
*/
void r3j000(int j2, int j3, C6 *c6);
-/*! \brief C++ porting of R3JJR
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients for JJ transitions.
+ *
+ * This function calculates the 3j(J,J2,J3;-M2-M3,M2,M3) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j2: `int`
* \param j3: `int`
@@ -200,7 +245,10 @@ void r3j000(int j2, int j3, C6 *c6);
*/
void r3jjr(int j2, int j3, int m2, int m3, C6 *c6);
-/*! \brief C++ porting of R3JMR
+/*! \brief Compute the 3j symbol for Clebsch-Gordan coefficients for JM transitions.
+ *
+ * This function calculates the 3j(J,J2,J3;M1,M,-M1-M) symbol for the Clebsch-Gordan
+ * coefficients. See Appendix a.3.1 in Borghese, Denti & Saija (2007).
*
* \param j1: `int`
* \param j2: `int`
@@ -210,7 +258,12 @@ void r3jjr(int j2, int j3, int m2, int m3, C6 *c6);
*/
void r3jmr(int j1, int j2, int j3, int m1, C6 *c6);
-/*! \brief C++ porting of RABA
+/*! \brief Compute radiation torques on a particle in Cartesian coordinates.
+ *
+ * This function computes radiation torque on on a cluster of spheres as the
+ * result of the difference between the extinction and the scattering
+ * contributions for a Cartesian coordinate system, as `rabas()`. See Sec. 4.9
+ * in Borghese, Denti & Saija (2007).
*
* \param le: `int`
* \param am0m: Matrix of complex.
@@ -225,7 +278,10 @@ void raba(
std::complex<double> **tqcpe, double **tqcs, std::complex<double> **tqcps
);
-/*! \brief C++ porting of RFTR
+/*! \brief Compute the radiation force Cartesian components.
+ *
+ * This function computes the Cartesian components of the radiation force
+ * exerted on a particle. See Sec. 3.2.1 in Borghese, Denti & Saija (2007).
*
* \param u: `double *`
* \param up: `double *`
@@ -248,7 +304,12 @@ void rftr(
double &fy, double &fz
);
-/*! \brief C++ porting of SCR0
+/*! \brief Compute Mie cross-sections for the sphere units in the cluster.
+ *
+ * This function computes the scattering, absorption and extinction cross-sections
+ * for the spheres composing the cluster, in terms of Mie coefficients, together
+ * with the Forward Scattering Amplitude. See Sec. 4.2.1 in Borghese, Denti & Saija
+ * (2007).
*
* \param vk: `double` Wave number
* \param exri: `double` External medium refractive index.
@@ -259,7 +320,10 @@ void rftr(
*/
void scr0(double vk, double exri, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 * c4);
-/*! \brief C++ porting of SCR2
+/*! \brief Compute the scattering amplitude for a single sphere in an aggregate.
+ *
+ * This function computes the scattering amplitude for single spheres constituting
+ * an aggregate. See Sec. 4.2.1 in Borghese, Denti & Saija (2007).
*
* \param vk: `double` Wave number.
* \param vkarg: `double` QUESTION: definition?
@@ -275,7 +339,11 @@ void scr2(
C3 *c3, C4 *c4
);
-/*! \brief C++ porting of STR
+/*! \brief Transform sphere Cartesian coordinates to spherical coordinates.
+ *
+ * This function transforms the Cartesian coordinates of the spheres in an aggregate
+ * to radial coordinates, then it calls `sphar()` to calculate the vector of spherical
+ * harmonics of the incident field.
*
* \param rcf: `double **` Matrix of double.
* \param c1: `C1 *` Pointer to a C1 instance.
@@ -286,7 +354,12 @@ void scr2(
*/
void str(double **rcf, C1 *c1, C1_AddOns *c1ao, C3 *c3, C4 *c4, C6 *c6);
-/*! \brief C++ porting of TQR
+/*! \brief Compute radiation torques on particles on a k-vector oriented system.
+ *
+ * This function computes the radiation torques resulting from the difference
+ * between absorption and scattering contributions, like `rabas()`, but for a
+ * coordinate system oriented along the wave vector and its orthogonal axis. See
+ * Sec. 4.9 in Borghese, Denti & Saija (2007).
*
* \param u: `double *`
* \param up: `double *`
@@ -305,7 +378,10 @@ void tqr(
double &ten, double &tek, double &tsp, double &tsn, double &tsk
);
-/*! \brief C++ porting of ZTM
+/*! \brief Calculate the single-centered inversion of the M-matrix.
+ *
+ * This function computes the single-centered inverrted M-matrix appearing in Eq. (5.28)
+ * of Borghese, Denti & Saija (2007).
*
* \param am: Matrix of complex.
* \param c1: `C1 *` Pointer to a C1 instance.
diff --git a/src/include/sph_subs.h b/src/include/sph_subs.h
index 9526b34341f85fc4a8be09a3ba259b655a05982a..2aafef2233caa43a8f76377afc65db380f056c45 100644
--- a/src/include/sph_subs.h
+++ b/src/include/sph_subs.h
@@ -14,7 +14,7 @@
#ifndef INCLUDE_SPH_SUBS_H_
#define INCLUDE_SPH_SUBS_H_
-/*! \brief Compute the asymmetry-corrected scattering cross-section.
+/*! \brief Compute the asymmetry-corrected scattering cross-section of a sphere.
*
* This function computes the product between the geometrical asymmetry parameter and
* the scattering cross-section. See Sec. 3.2.1 of Borghese, Denti & Saija (2007).
@@ -179,9 +179,9 @@ void pwma(
int isq, C1 *c1
);
-/*! \brief Compute radiation torques on particles.
+/*! \brief Compute radiation torques on a single sphere in Cartesian coordinates.
*
- * This function computes radiation torque on the particle as a result
+ * This function computes radiation torque on a sphere unit as the result
* of the difference between the extinction and the scattering contributions.
* See Sec. 4.9 in Borghese, Denti & Saija (2007).
*
@@ -321,7 +321,12 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1);
*/
void thdps(int lm, double ****zpv);
-/*! \brief C++ porting of UPVMP
+/*! \brief Compute the unitary vectors onb the polarization plane and its orthogonal
+ * direction.
+ *
+ * This function computes the unitary vectors lying on the polarization plane and on
+ * its orthogonal direction, to optimize the identification of the scattering geometry.
+ * See Sec. 2.3 in Borghese, Denti & Saija (2007).
*
* \param thd: `double`
* \param phd: `double`
diff --git a/src/libnptm/tra_subs.cpp b/src/libnptm/tra_subs.cpp
index c9a490c446b7951ac1072c810506395841d46c42..8aa135255e1968a8d3141be9cefe1cdd4cf6f797 100644
--- a/src/libnptm/tra_subs.cpp
+++ b/src/libnptm/tra_subs.cpp
@@ -362,7 +362,7 @@ void sampoa(
}
void wamff(
- complex<double> *wk, double x, double y, double z, int lm, double apfafa,
+ std::complex<double> *wk, double x, double y, double z, int lm, double apfafa,
double tra, double spd, double rir, double ftcn, int lmode, double pmf
) {
const int nlmm = lm * (lm + 2);