/*! \file sph_subs.h
 *
 * \brief C++ porting of SPH functions and subroutines.
 *
 * This library includes a collection of functions that are used to solve the
 * scattering problem in the case of a single sphere. Some of these functions
 * are also generalized to the case of clusters of spheres. In most cases, the
 * results of calculations do not fall back to fundamental data types. They are
 * rather multi-component structures. In order to manage access to such variety
 * of return values, most functions are declared as `void` and they operate on
 * output arguments passed by reference.
 */

#ifndef INCLUDE_SPH_SUBS_H_
#define INCLUDE_SPH_SUBS_H_

/*! \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 ****` Geometrical asymmetry parameter coefficients matrix.
 * \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 *` Geometrical asymmetry parameter-corrected cross-section.
 */
void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps);

/*! \brief Complex Bessel Function.
 *
 * This function computes the complex spherical Bessel funtions \f$j\f$. It uses the
 * auxiliary functions `msta1()` and `msta2()` to determine the starting point for
 * backward recurrence. This is the `CSPHJ` implementation of the `specfun` library.
 *
 * \param n: `int` Order of the function.
 * \param z: `complex<double>` Argumento of the function.
 * \param nm: `int &` Highest computed order.
 * \param csj: Vector of complex. The desired function \f$j\f$.
 */
void cbf(int n, std::complex<double> z, int &nm, std::complex<double> csj[]);

/*! \brief Clebsch-Gordan coefficients.
 *
 * This function comutes the Clebsch-Gordan coefficients for the irreducible spherical
 * tensors. See Sec. 1.6.3 and Table 1.1 in Borghese, Denti & Saija (2007).
 *
 * \param lmpml: `int`
 * \param mu: `int`
 * \param l: `int`
 * \param m: `int`
 * \return result: `double` Clebsh-Gordan coefficient.
 */
double cg1(int lmpml, int mu, int l, int m);

/*! \brief Conjugate of a double precision complex number
 *
 * \param z: `complex<double>` The input complex number.
 * \return result: `complex<double>` The conjugate of the input number.
 */
std::complex<double> dconjg(std::complex<double> z);

/*! \brief Comute the continuous variation of the refractive index and of its radial derivative.
 *
 * This function implements the continuous variation of the refractive index and of its radial
 * derivative through the materials that constitute the sphere and its surrounding medium. See
 * Sec. 5.2 in Borghese, Denti & Saija (2007).
 *
 * \param npntmo: `int`
 * \param ns: `int`
 * \param i: `int`
 * \param ic: `int`
 * \param vk: `double`
 * \param c1: `C1 *` Pointer to `C1` data structure.
 * \param c2: `C2 *` Pointer to `C2` data structure.
 */
void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2);

/*! \brief Compute Mie scattering coefficients.
 *
 * This function determines the L-dependent Mie scattering coefficients \f$a_l\f$ and \f$b_l\f$
 * for the cases of homogeneous spheres, radially non-homogeneous spheres and, in case of sphere
 * with dielectric function, sustaining longitudinal waves. See Sec. 5.1 in Borghese, Denti
 * & Saija (2007).
 *
 * \param li: `int` Maximum field expansion order.
 * \param i: `int`
 * \param npnt: `int`
 * \param npntts: `int`
 * \param vk: `double` Wave number in scale units.
 * \param exdc: `double` External medium dielectric constant.
 * \param exri: `double` External medium refractive index.
 * \param c1: `C1 *` Pointer to a `C1` data structure.
 * \param c2: `C2 *` Pointer to a `C2` data structure.
 * \param jer: `int &` Reference to integer error code variable.
 * \param lcalc: `int &` Reference to integer variable recording the maximum expansion order accounted for.
 * \param arg: `complex<double> &`
 */
void dme(
	 int li, int i, int npnt, int npntts, double vk, double exdc, double exri,
	 C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg
);

/*! \brief Bessel function calculation control parameters.
 *
 * This function determines the control parameters required to calculate the Bessel
 * functions up to the desired precision.
 *
 * \param n: `int`
 * \param x: `double`
 * \return result: `double`
 */
double envj(int n, double x);

/*! \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 **`
 * \param cmul: `double **`
 */
void mmulc(std::complex<double> *vint, double **cmullr, double **cmul);

/*! \brief Starting point for Bessel function magnitude.
 *
 * This function determines the starting point for backward recurrence such that
 * the magnitude of all \f$J\f$ and \f$j\f$ functions is of the order of \f$10^{-mp}\f$.
 *
 * \param x: `double` Absolute value of the argumetn to \f$J\f$ or \f$j\f$.
 * \param mp: `int` Requested order of magnitude.
 * \return result: `int` The necessary starting point.
 */
int msta1(double x, int mp);

/*! \brief Starting point for Bessel function precision.
 *
 * This function determines the starting point for backward recurrence such that
 * all \f$J\f$ and \f$j\f$ functions have `mp` significant digits.
 *
 * \param x: `double` Absolute value of the argumetn to \f$J\f$ or \$j\f$.
 * \param n: `int` Order of the function.
 * \param mp: `int` Requested number of significant digits.
 * \return result: `int` The necessary starting point.
 */
int msta2(double x, int n, int mp);

/*! \brief Compute the amplitude of the orthogonal unit vector.
 *
 * This function computes the amplitude of the orthogonal unit vector for a geometry
 * based either on the scattering plane or on the meridional plane. It is used by `upvsp()`
 * and `upvmp()`. See Sec. 2.7 in Borghese, Denti & Saija (2007).
 *
 * \param u1: `double *`
 * \param u2: `double *`
 * \param u3: `double *`
 * \param iorth: `int`
 * \param torth: `double`
 */
void orunve(double *u1, double *u2, double *u3, int iorth, double torth);

/*! \brief Compute incident and scattered field amplitudes.
 *
 * This function computes the amplitudes of the incident and scattered field on the
 * basis of the multi-polar expansion. See Sec. 4.1.1 in Borghese, Denti and Saija (2007).
 *
 * \param up: `double *`
 * \param un: `double *`
 * \param ylm: Vector of complex. Field polar spherical harmonics.
 * \param inpol: `int` Incident field polarization type (0 - linear; 1 - circular).
 * \param lw: `int`
 * \param isq: `int`
 * \param c1: `C1 *` Pointer to a `C1` data structure.
 */
void pwma(
	  double *up, double *un, std::complex<double> *ylm, int inpol, int lw,
	  int isq, C1 *c1
);

/*! \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` 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.
 * \param tqsps: Matrix of complex.
 */
void rabas(
	   int inpol, int li, int nsph, C1 *c1, double **tqse, std::complex<double> **tqspe,
	   double **tqss, std::complex<double> **tqsps
);

/*! \brief Real Bessel Function.
 *
 * This function computes the real spherical Bessel funtions \f$j\f$. It uses the
 * auxiliary functions `msta1()` and `msta2()` to determine the starting point for
 * backward recurrence. This is the `SPHJ` implementation of the `specfun` library.
 *
 * \param n: `int` Order of the function.
 * \param x: `double` Argument of the function.
 * \param nm: `int &` Highest computed order.
 * \param sj: `double[]` The desired function \f$j\f$.
 */
void rbf(int n, double x, int &nm, double sj[]);

/*! \brief Soft layer radial function and derivative.
 *
 * This function determines the radial function and its derivative for a soft layer
 * in a radially non homogeneous sphere. See Sec. 5.1 in Borghese, Denti & Saija (2007).
 *
 * \param npntmo: `int`
 * \param step: `double`
 * \param dcc: `complex<double>`
 * \param x: `double &`
 * \param lpo: `int`
 * \param y1: `complex<double> &`
 * \param y2: `complex<double> &`
 * \param dy1: `complex<double> &`
 * \param dy2: `complex<double> &`
 */
void rkc(
	 int npntmo, double step, std::complex<double> dcc, double &x, int lpo,
	 std::complex<double> &y1, std::complex<double> &y2, std::complex<double> &dy1,
	 std::complex<double> &dy2
);

/*! \brief Transition layer radial function and derivative.
 *
 * This function determines the radial function and its derivative for a transition layer
 * in a radially non homogeneous sphere. See Sec. 5.1 in Borghese, Denti & Saija (2007).
 *
 * \param npntmo: `int`
 * \param step: `double`
 * \param x: `double &`
 * \param lpo: `int`
 * \param y1: `complex<double> &`
 * \param y2: `complex<double> &`
 * \param dy1: `complex<double> &`
 * \param dy2: `complex<double> &`
 * \param c2: `C2 *` Pointer to a `C2` data structure.
 */
void rkt(
	 int npntmo, double step, double &x, int lpo, std::complex<double> &y1,
	 std::complex<double> &y2, std::complex<double> &dy1, std::complex<double> &dy2,
	 C2 *c2
);

/*! \brief Spherical Bessel functions.
 *
 * This function computes the spherical Bessel functions \f$y\f$. It adopts the `SPHJY`
 * implementation of the `specfun` library.
 *
 * \param n: `int` Order of the function (from 0 up).
 * \param x: `double` Argumento of the function (\f$x > 0\f$).
 * \param nm: `int &` Highest computed order.
 * \param sy: `double[]` The desired function \f$y\f$.
 */
void rnf(int n, double x, int &nm, double sy[]);

/*! \brief Spherical harmonics for given direction.
 *
 * This function computes the field spherical harmonics for a given direction. See Sec.
 * 1.5.2 in Borghese, Denti & Saija (2007).
 *
 * \param cosrth: `double` Cosine of direction's elevation.
 * \param sinrth: `double` Sine of direction's elevation.
 * \param cosrph: `double` Cosine of direction's azimuth.
 * \param sinrph: `double` Sine of direction's azimuth.
 * \param ll: `int` L value expansion order.
 * \param ylm: Vector of complex. The requested spherical harmonics.
 */
void sphar(
	   double cosrth, double sinrth, double cosrph, double sinrph,
	   int ll, std::complex<double> *ylm
);

/*! \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` 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);

/*! \brief C++ Compute the scattering amplitude and the scattered field intensity.
 *
 * 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);

/*! \brief Determine the geometrical asymmetry parameter coefficients.
 *
 * 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);

/*! \brief C++ porting of UPVMP
 *
 * \param thd: `double`
 * \param phd: `double`
 * \param icspnv: `int`
 * \param cost: `double`
 * \param sint: `double`
 * \param cosp: `double`
 * \param sinp: `double`
 * \param u: `double *`
 * \param up: `double *`
 * \param un: `double *`
 */
void upvmp(
	   double thd, double phd, int icspnv, double &cost, double &sint,
	   double &cosp, double &sinp, double *u, double *up, double *un
);

/*! \brief Compute the unitary vector perpendicular to incident and scattering plane.
 *
 * This function computes the unitary vector perpendicular to the incident and scattering
 * plane in a geometry based on the scattering plane. It uses `orunve()`. See Sec. 2.7 in
 * Borghese, Denti & Saija (2007).
 *
 * \param u: `double *`
 * \param upmp: `double *`
 * \param unmp: `double *`
 * \param us: `double *`
 * \param upsmp: `double *`
 * \param unsmp: `double *`
 * \param up: `double *`
 * \param un: `double *`
 * \param ups: `double *`
 * \param uns: `double *`
 * \param duk: `double *`
 * \param isq: `int &`
 * \param ibf: `int &`
 * \param scand: `double &` Scattering angle in degrees.
 * \param cfmp: `double &`
 * \param sfmp: `double &`
 * \param cfsp: `double &`
 * \param sfsp: `double &`
 */
void upvsp(
	   double *u, double *upmp, double *unmp, double *us, double *upsmp, double *unsmp,
	   double *up, double *un, double *ups, double *uns, double *duk, int &isq,
	   int &ibf, double &scand, double &cfmp, double &sfmp, double &cfsp, double &sfsp
);

/*! \brief Compute meridional plane-referred geometrical asymmetry parameter coefficients.
 *
 * This function computes the coeffcients that define the geometrical symmetry parameter
 * as defined with respect to the meridional plane. It makes use of `sphar()` and `pwma()`.
 * See Sec. 3.2.1 in Borghese, Denti & Saija (2007).
 *
 * \param iis: `int`
 * \param cost: `double` Cosine of the elevation angle.
 * \param sint: `double` Sine of the elevation angle.
 * \param cosp: `double` Cosine of the azimuth angle.
 * \param sinp: `double` Sine of the azimuth angle.
 * \param inpol: `int` Incident field polarization type (0 - linear; 1 - circular).
 * \param lm: `int` Maximum field expansion orde.
 * \param idot: `int`
 * \param nsph: `int` Number of spheres.
 * \param arg: `double *`
 * \param u: `double *`
 * \param up: `double *`
 * \param un: `double *`
 * \param c1: `C1 *` Pointer to a `C1` data structure.
 */
void wmamp(
	   int iis, double cost, double sint, double cosp, double sinp, int inpol,
	   int lm, int idot, int nsph, double *arg, double *u, double *up,
	   double *un, C1 *c1
);

/*! \brief Compute the scattering plane-referred geometrical asymmetry parameter coefficients.
 *
 * This function computes the coefficients that define the geometrical asymmetry parameter based
 * on the L-value with respect to the scattering plane. It uses `sphar()` and `pwma()`. See Sec.
 * 3.2.1 in Borghese, Denti and Saija (2007).
 *
 * \param cost: `double` Cosine of elevation angle.
 * \param sint: `double` Sine of elevation angle.
 * \param cosp: `double` Cosine of azimuth angle.
 * \param sinp: `double` Sine of azimuth angle.
 * \param costs: `double` Cosine of scattering elevation angle.
 * \param sints: `double` Sine of scattering elevation angle.
 * \param cosps: `double` Cosine of scattering azimuth angle.
 * \param sinps: `double` Sine of scattering azimuth angle.
 * \param u: `double *`
 * \param up: `double *`
 * \param un: `double *`
 * \param us: `double *`
 * \param ups: `double *`
 * \param uns: `double *`
 * \param isq: `int`
 * \param ibf: `int`
 * \param inpol: `int` Incident field polarization (0 - linear; 1 -circular).
 * \param lm: `int` Maximum field expansion order.
 * \param idot: `int`
 * \param nsph: `int` Number opf spheres.
 * \param argi: `double *`
 * \param args: `double *`
 * \param c1: `C1 *` Pointer to a `C1` data structure.
 */
void wmasp(
	   double cost, double sint, double cosp, double sinp, double costs, double sints,
	   double cosps, double sinps, double *u, double *up, double *un, double *us,
	   double *ups, double *uns, int isq, int ibf, int inpol, int lm, int idot,
	   int nsph, double *argi, double *args, C1 *c1
);

#endif /* SRC_INCLUDE_SPH_SUBS_H_ */
