Add documentation for CLUSTER routines

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.

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`

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);