Expand documentation of SPH subroutines

src/include/sph_subs.h

@@ -96,14 +96,17 @@ double cg1(int lmpml, int mu, int l, int m) {
   return result;
 }
 
-/*! \brief C++ porting of APS
+/*! \brief Compute the asymmetry-corrected scattering cross-section.
+ *
+ * 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).
  *
  * \param zpv: `double ****`
- * \param li: `int`
- * \param nsph: `int`
- * \param c1: `C1 *`
+ * \param li: `int` Maximum field expansion order.
+ * \param nsph: `int` Number of spheres.
+ * \param c1: `C1 *` Pointer to a `C1` data structure.
  * \param sqk: `double`
- * \param gaps: `double *`
+ * \param gaps: `double *` Geometrical asymmetry parameter-corrected cross-section.
  */
 void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) {
   std::complex<double> cc0 = std::complex<double>(0.0, 0.0);
@@ -338,7 +341,10 @@ void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]) {
   }
 }
 
-/*! \brief C++ porting of MMULC
+/*! \brief Compute the Mueller Transformation Matrix.
+ *
+ * This function computes the Mueller Transformation Matrix, or Phase Matrix. See
+ * Sec. 2.8.1 of Borghese, Denti & Saija (2007).
  *
  * \param vint: Vector of complex.
  * \param cmullr: `double **`
@@ -510,15 +516,16 @@ void pwma(
   }
 }
 
-/*! \brief C++ porting of RABAS
+/*! \brief Compute radiation torques on particles.
  *
  * This function computes radiation torque on the particle as a result
  * of the difference between the extinction and the scattering contributions.
+ * See Sec. 4.9 in Borghese, Denti & Saija (2007).
  *
- * \param inpol: `int`
- * \param li: `int`
- * \param nsph: `int`
- * \param c1: `C1 *`
+ * \param inpol: `int` Incident polarization type (0 - linear; 1 - circular)
+ * \param li: `int` Maximum field expansion order.
+ * \param nsph: `int` Number of spheres.
+ * \param c1: `C1 *` Pointer to `C1` data structure.
  * \param tqse: Matrix of complex.
  * \param tqspe: Matrix of complex.
  * \param tqss: Matrix of complex.
@@ -809,14 +816,17 @@ void rnf(int n, double x, int &nm, double sy[]) {
   return;
 }
 
-/*! \brief C++ porting of SSCR0
+/*! \brief Compute scattering, absorption and extinction cross-sections.
+ *
+ * This function computes the scattering, absorption and extinction cross-sections in terms
+ * of Forward Scattering Amplitudes. See Sec. 4.2.1 in Borghese, Denti & Saija (2007).
  *
  * \param tfsas: `complex<double> &`
- * \param nsph: `int`
- * \param lm: `int`
- * \param vk: `double`
- * \param exri: `double`
- * \param c1: `C1 *`
+ * \param nsph: `int` Number of spheres.
+ * \param lm: `int` Maximum field expansion order.
+ * \param vk: `double` Wave number in scale units.
+ * \param exri: `double` External medium refractive index.
+ * \param c1: `C1 *` Pointer to a `C1` data structure.
  */
 void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri, C1 *c1) {
   std::complex<double> sum21, rm, re, csam;
@@ -859,13 +869,16 @@ void sscr0(std::complex<double> &tfsas, int nsph, int lm, double vk, double exri
   }
 }
 
-/*! \brief C++ porting of SSCR2
+/*! \brief C++ Compute the scattering amplitude and the scattered field intensity.
  *
- * \param nsph: `int`
- * \param lm: `int`
- * \param vk: `double`
- * \param exri: `double`
- * \param c1: `C1 *`
+ * The role of this function is to compute the scattering amplitude and the intensity of
+ * the scattered field. See Sec. 4.2 in Borghese, Denti & Saija (2007).
+ *
+ * \param nsph: `int` Number of spheres.
+ * \param lm: `int` Maximum field expansion order.
+ * \param vk: `double` Wave number in scale units.
+ * \param exri: `double` External medium refractive index.
+ * \param c1: `C1 *` Pointer to a `C1` data structure.
  */
 void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) {
   std::complex<double> s11, s21, s12, s22, rm, re, csam, cc;
@@ -1011,10 +1024,14 @@ void sphar(
   goto label40;
 }
 
-/*! \brief C++ porting of THDPS
+/*! \brief Determine the geometrical asymmetry parameter coefficients.
  *
- * \param lm: `int`
- * \param zpv: `double ****`
+ * This function computes the coefficients that enter the definition of the geometrical
+ * asymmetry parameter based on the L-value of the field expansion order. See Sec. 3.2.1
+ * in Borghese, Denti & Saija (2007).
+ *
+ * \param lm: `int` Maximum field expansion order.
+ * \param zpv: `double ****` Matrix of geometrical asymmetry parameter coefficients.
  */
 void thdps(int lm, double ****zpv) {
   //for (int l10 = 0; l10 < lm; l10++) { // 0-init, can be omitted