Loading src/include/sph_subs.h +157 −105 Original line number Diff line number Diff line Loading @@ -2,13 +2,13 @@ * * \brief C++ porting of SPH functions and subroutines. * * Remember that FORTRAN passes arguments by reference, so, every time we use * a subroutine call, we need to add a referencing layer to the C++ variable. * All the functions defined below need to be properly documented and ported * to C++. * * Currently, only basic documenting information about functions and parameter * types are given, to avoid doxygen warning messages. * This library includes a collection of function 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_COMMONS_H_ Loading @@ -22,8 +22,8 @@ /*! \brief Conjugate of a double precision complex number * * \param z: `std::complex\<double\>` The input complex number. * \return result: `std::complex\<double\>` The conjugate of the input 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) { double zreal = z.real(); Loading @@ -31,13 +31,16 @@ std::complex<double> dconjg(std::complex<double> z) { return std::complex<double>(zreal, -zimag); } /*! \brief C++ porting of CG1 /*! \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` * \return result: `double` Clebsh-Gordan coefficient. */ double cg1(int lmpml, int mu, int l, int m) { double result = 0.0; Loading Loading @@ -101,7 +104,7 @@ double cg1(int lmpml, int mu, int l, int m) { * 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 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. Loading Loading @@ -158,15 +161,19 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) { } } /*! \brief C++ porting of DIEL /*! \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 *` * \param c2: `C2 *` * \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) { const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1]; Loading @@ -187,7 +194,10 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) { } } /*! \brief C++ porting of ENVJ /*! \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` Loading @@ -205,11 +215,14 @@ double envj(int n, double x) { return result; } /*! \brief C++ porting of MSTA1 /*! \brief Starting point for Bessel function magnitude. * * \param x: `double` * \param mp: `int` * \return result: `int` * This function determines the starting point for backward recurrence such that * the magnitude of all \f$J\f$ and \$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 \$j\f$. * \param mp: `int` Requested order of magnitude. * \return result: `int` The necessary starting point. */ int msta1(double x, int mp) { int result = 0; Loading @@ -236,12 +249,15 @@ int msta1(double x, int mp) { return result; } /*! \brief C++ porting of MSTA2 /*! \brief Starting point for Bessel function precision. * * \param x: `double` * \param n: `int` * \param mp: `int` * \return result: `int` * This function determines the starting point for backward recurrence such that * all \f$J\f$ and \$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) { int result = 0; Loading Loading @@ -276,14 +292,16 @@ int msta2(double x, int n, int mp) { return result; } /*! \brief C++ porting of CBF /*! \brief Complex Bessel Function. * * This is the 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` * \param z: `complex\<double\>` * \param nm: `int &` * \param csj: Vector of complex. * \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[]) { /* Loading Loading @@ -401,7 +419,11 @@ void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) { cmul[3][3] = s21 - s34; } /*! \brief C++ porting of ORUNVE /*! \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 *` Loading @@ -427,15 +449,18 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) { u3[2] = u1[0] * u2[1] - u1[1] * u2[0]; } /*! \brief C++ porting of PWMA /*! \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 * \param inpol: `int` * \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 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void pwma( double *up, double *un, std::complex<double> *ylm, int inpol, int lw, Loading Loading @@ -587,14 +612,16 @@ void rabas( } } /*! \brief C++ porting of RBF /*! \brief Real Bessel Function. * * This is the 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` * \param x: `double` * \param nm: `int &` * \param sj: `double[]` * \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[]) { /* Loading Loading @@ -652,17 +679,20 @@ void rbf(int n, double x, int &nm, double sj[]) { } } /*! \brief C++ porting of RKC /*! \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 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\> &` * \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, Loading Loading @@ -702,17 +732,20 @@ void rkc( } } /*! \brief C++ porting of RKT /*! \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 *` * \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, Loading Loading @@ -763,14 +796,15 @@ void rkt( } } /*! \brief C++ porting of RNF /*! \brief Spherical Bessel functions. * * This is a real spherical Bessel function. * This function computes the spherical Bessel functions \f$y\$. It adopts the `SPHJY` * implementation of the `specfun` library. * * \param n: `int` * \param x: `double` * \param nm: `int &` * \param sy: `double[]` * \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[]) { /* Loading Loading @@ -936,14 +970,17 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) { } } /*! \brief C++ porting of SPHAR /*! \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` * \param sinrth: `double` * \param cosrph: `double` * \param sinrph: `double` * \param ll: `int` * \param ylm: Vector of complex. * \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, Loading Loading @@ -1113,7 +1150,11 @@ void upvmp( } } /*! \brief C++ porting of UPVSP /*! \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 *` Loading @@ -1128,7 +1169,7 @@ void upvmp( * \param duk: `double *` * \param isq: `int &` * \param ibf: `int &` * \param scand: `double &` * \param scand: `double &` Scattering angle in degrees. * \param cfmp: `double &` * \param sfmp: `double &` * \param cfsp: `double &` Loading Loading @@ -1202,22 +1243,26 @@ void upvsp( sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2]; } /*! \brief C++ porting of WMAMP /*! \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` * \param sint: `double` * \param cosp: `double` * \param sinp: `double` * \param inpol: `int` * \param lm: `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` * \param nsph: `int` Number of spheres. * \param arg: `double *` * \param u: `double *` * \param up: `double *` * \param un: `double *` * \param c1: `C1 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void wmamp( int iis, double cost, double sint, double cosp, double sinp, int inpol, Loading @@ -1243,21 +1288,24 @@ void wmamp( } } sphar(cost, sint, cosp, sinp, lm, ylm); //printf("DEBUG: in WMAMP and calling PWMA with lm = %d and iis = %d\n", lm, iis); pwma(up, un, ylm, inpol, lm, iis, c1); delete[] ylm; } /*! \brief C++ porting of WMASP /*! \brief Compute the scattering plane-referred geometrical asymmetry parameter coefficients. * * \param cost: `double` * \param sint: `double` * \param cosp: `double` * \param sinp: `double` * \param costs: `double` * \param sints: `double` * \param cosps: `double` * \param sinps: `double` * 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 *` Loading @@ -1266,13 +1314,13 @@ void wmamp( * \param uns: `double *` * \param isq: `int` * \param ibf: `int` * \param inpol: `int` * \param lm: `int` * \param inpol: `int` Incident field polarization (0 - linear; 1 -circular). * \param lm: `int` Maximum field expansion order. * \param idot: `int` * \param nsph: `int` * \param nsph: `int` Number opf spheres. * \param argi: `double *` * \param args: `double *` * \param c1: `C1 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void wmasp( double cost, double sint, double cosp, double sinp, double costs, double sints, Loading Loading @@ -1306,34 +1354,38 @@ void wmasp( } } sphar(cost, sint, cosp, sinp, lm, ylm); //printf("DEBUG: in WMASP and calling PWMA with isq = %d\n", isq); pwma(up, un, ylm, inpol, lm, isq, c1); if (ibf >= 0) { sphar(costs, sints, cosps, sinps, lm, ylm); //printf("DEBUG: in WMASP and calling PWMA with isq = 2 and ibf = %d\n", ibf); pwma(ups, uns, ylm, inpol, lm, 2, c1); } delete[] ylm; } /*! \brief C++ porting of DME /*! \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` * \param li: `int` Maximum field expansion order. * \param i: `int` * \param npnt: `int` * \param npntts: `int` * \param vk: `double` * \param exdc: `double` * \param exri: `double` * \param c1: `C1 *` * \param c2: `C2 *` * \param jer: `int &` * \param lcalc: `int &` * \param arg: `complex<double> &`. * \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) { C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg ) { const int lipo = li + 1; const int lipt = li + 2; double *rfj = new double[lipt]; Loading src/include/tra_subs.h +23 −3 Original line number Diff line number Diff line Loading @@ -7,12 +7,32 @@ #endif // Structures for TRAPPING (that will probably move to Commons) /*! \brief CIL data structure. * * A structure containing field expansion order configuration. */ struct CIL { int le, nlem, nlemt, mxmpo, mxim; //! Maximum L expansion of the electric field. int le; //! le * (le + 1). int nlem; //! 2 * nlem. int nlemt; //! Maximum field expansion order + 1. int mxmpo; //! 2 * mxmpo - 1. int mxim; }; /*! \brief CCR data structure. * * A structure containing geometrical asymmetry parameter normalization coefficients. */ struct CCR { double cof, cimu; //! First coefficient. double cof; //! Second coefficient. double cimu; }; //End of TRAPPING structures Loading Loading
src/include/sph_subs.h +157 −105 Original line number Diff line number Diff line Loading @@ -2,13 +2,13 @@ * * \brief C++ porting of SPH functions and subroutines. * * Remember that FORTRAN passes arguments by reference, so, every time we use * a subroutine call, we need to add a referencing layer to the C++ variable. * All the functions defined below need to be properly documented and ported * to C++. * * Currently, only basic documenting information about functions and parameter * types are given, to avoid doxygen warning messages. * This library includes a collection of function 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_COMMONS_H_ Loading @@ -22,8 +22,8 @@ /*! \brief Conjugate of a double precision complex number * * \param z: `std::complex\<double\>` The input complex number. * \return result: `std::complex\<double\>` The conjugate of the input 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) { double zreal = z.real(); Loading @@ -31,13 +31,16 @@ std::complex<double> dconjg(std::complex<double> z) { return std::complex<double>(zreal, -zimag); } /*! \brief C++ porting of CG1 /*! \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` * \return result: `double` Clebsh-Gordan coefficient. */ double cg1(int lmpml, int mu, int l, int m) { double result = 0.0; Loading Loading @@ -101,7 +104,7 @@ double cg1(int lmpml, int mu, int l, int m) { * 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 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. Loading Loading @@ -158,15 +161,19 @@ void aps(double ****zpv, int li, int nsph, C1 *c1, double sqk, double *gaps) { } } /*! \brief C++ porting of DIEL /*! \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 *` * \param c2: `C2 *` * \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) { const double dif = c1->rc[i - 1][ns] - c1->rc[i - 1][ns - 1]; Loading @@ -187,7 +194,10 @@ void diel(int npntmo, int ns, int i, int ic, double vk, C1 *c1, C2 *c2) { } } /*! \brief C++ porting of ENVJ /*! \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` Loading @@ -205,11 +215,14 @@ double envj(int n, double x) { return result; } /*! \brief C++ porting of MSTA1 /*! \brief Starting point for Bessel function magnitude. * * \param x: `double` * \param mp: `int` * \return result: `int` * This function determines the starting point for backward recurrence such that * the magnitude of all \f$J\f$ and \$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 \$j\f$. * \param mp: `int` Requested order of magnitude. * \return result: `int` The necessary starting point. */ int msta1(double x, int mp) { int result = 0; Loading @@ -236,12 +249,15 @@ int msta1(double x, int mp) { return result; } /*! \brief C++ porting of MSTA2 /*! \brief Starting point for Bessel function precision. * * \param x: `double` * \param n: `int` * \param mp: `int` * \return result: `int` * This function determines the starting point for backward recurrence such that * all \f$J\f$ and \$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) { int result = 0; Loading Loading @@ -276,14 +292,16 @@ int msta2(double x, int n, int mp) { return result; } /*! \brief C++ porting of CBF /*! \brief Complex Bessel Function. * * This is the 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` * \param z: `complex\<double\>` * \param nm: `int &` * \param csj: Vector of complex. * \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[]) { /* Loading Loading @@ -401,7 +419,11 @@ void mmulc(std::complex<double> *vint, double **cmullr, double **cmul) { cmul[3][3] = s21 - s34; } /*! \brief C++ porting of ORUNVE /*! \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 *` Loading @@ -427,15 +449,18 @@ void orunve( double *u1, double *u2, double *u3, int iorth, double torth) { u3[2] = u1[0] * u2[1] - u1[1] * u2[0]; } /*! \brief C++ porting of PWMA /*! \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 * \param inpol: `int` * \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 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void pwma( double *up, double *un, std::complex<double> *ylm, int inpol, int lw, Loading Loading @@ -587,14 +612,16 @@ void rabas( } } /*! \brief C++ porting of RBF /*! \brief Real Bessel Function. * * This is the 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` * \param x: `double` * \param nm: `int &` * \param sj: `double[]` * \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[]) { /* Loading Loading @@ -652,17 +679,20 @@ void rbf(int n, double x, int &nm, double sj[]) { } } /*! \brief C++ porting of RKC /*! \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 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\> &` * \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, Loading Loading @@ -702,17 +732,20 @@ void rkc( } } /*! \brief C++ porting of RKT /*! \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 *` * \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, Loading Loading @@ -763,14 +796,15 @@ void rkt( } } /*! \brief C++ porting of RNF /*! \brief Spherical Bessel functions. * * This is a real spherical Bessel function. * This function computes the spherical Bessel functions \f$y\$. It adopts the `SPHJY` * implementation of the `specfun` library. * * \param n: `int` * \param x: `double` * \param nm: `int &` * \param sy: `double[]` * \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[]) { /* Loading Loading @@ -936,14 +970,17 @@ void sscr2(int nsph, int lm, double vk, double exri, C1 *c1) { } } /*! \brief C++ porting of SPHAR /*! \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` * \param sinrth: `double` * \param cosrph: `double` * \param sinrph: `double` * \param ll: `int` * \param ylm: Vector of complex. * \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, Loading Loading @@ -1113,7 +1150,11 @@ void upvmp( } } /*! \brief C++ porting of UPVSP /*! \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 *` Loading @@ -1128,7 +1169,7 @@ void upvmp( * \param duk: `double *` * \param isq: `int &` * \param ibf: `int &` * \param scand: `double &` * \param scand: `double &` Scattering angle in degrees. * \param cfmp: `double &` * \param sfmp: `double &` * \param cfsp: `double &` Loading Loading @@ -1202,22 +1243,26 @@ void upvsp( sfsp = uns[0] * upsmp[0] + uns[1] * upsmp[1] + uns[2] * upsmp[2]; } /*! \brief C++ porting of WMAMP /*! \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` * \param sint: `double` * \param cosp: `double` * \param sinp: `double` * \param inpol: `int` * \param lm: `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` * \param nsph: `int` Number of spheres. * \param arg: `double *` * \param u: `double *` * \param up: `double *` * \param un: `double *` * \param c1: `C1 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void wmamp( int iis, double cost, double sint, double cosp, double sinp, int inpol, Loading @@ -1243,21 +1288,24 @@ void wmamp( } } sphar(cost, sint, cosp, sinp, lm, ylm); //printf("DEBUG: in WMAMP and calling PWMA with lm = %d and iis = %d\n", lm, iis); pwma(up, un, ylm, inpol, lm, iis, c1); delete[] ylm; } /*! \brief C++ porting of WMASP /*! \brief Compute the scattering plane-referred geometrical asymmetry parameter coefficients. * * \param cost: `double` * \param sint: `double` * \param cosp: `double` * \param sinp: `double` * \param costs: `double` * \param sints: `double` * \param cosps: `double` * \param sinps: `double` * 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 *` Loading @@ -1266,13 +1314,13 @@ void wmamp( * \param uns: `double *` * \param isq: `int` * \param ibf: `int` * \param inpol: `int` * \param lm: `int` * \param inpol: `int` Incident field polarization (0 - linear; 1 -circular). * \param lm: `int` Maximum field expansion order. * \param idot: `int` * \param nsph: `int` * \param nsph: `int` Number opf spheres. * \param argi: `double *` * \param args: `double *` * \param c1: `C1 *` * \param c1: `C1 *` Pointer to a `C1` data structure. */ void wmasp( double cost, double sint, double cosp, double sinp, double costs, double sints, Loading Loading @@ -1306,34 +1354,38 @@ void wmasp( } } sphar(cost, sint, cosp, sinp, lm, ylm); //printf("DEBUG: in WMASP and calling PWMA with isq = %d\n", isq); pwma(up, un, ylm, inpol, lm, isq, c1); if (ibf >= 0) { sphar(costs, sints, cosps, sinps, lm, ylm); //printf("DEBUG: in WMASP and calling PWMA with isq = 2 and ibf = %d\n", ibf); pwma(ups, uns, ylm, inpol, lm, 2, c1); } delete[] ylm; } /*! \brief C++ porting of DME /*! \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` * \param li: `int` Maximum field expansion order. * \param i: `int` * \param npnt: `int` * \param npntts: `int` * \param vk: `double` * \param exdc: `double` * \param exri: `double` * \param c1: `C1 *` * \param c2: `C2 *` * \param jer: `int &` * \param lcalc: `int &` * \param arg: `complex<double> &`. * \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) { C1 *c1, C2 *c2, int &jer, int &lcalc, std::complex<double> &arg ) { const int lipo = li + 1; const int lipt = li + 2; double *rfj = new double[lipt]; Loading
src/include/tra_subs.h +23 −3 Original line number Diff line number Diff line Loading @@ -7,12 +7,32 @@ #endif // Structures for TRAPPING (that will probably move to Commons) /*! \brief CIL data structure. * * A structure containing field expansion order configuration. */ struct CIL { int le, nlem, nlemt, mxmpo, mxim; //! Maximum L expansion of the electric field. int le; //! le * (le + 1). int nlem; //! 2 * nlem. int nlemt; //! Maximum field expansion order + 1. int mxmpo; //! 2 * mxmpo - 1. int mxim; }; /*! \brief CCR data structure. * * A structure containing geometrical asymmetry parameter normalization coefficients. */ struct CCR { double cof, cimu; //! First coefficient. double cof; //! Second coefficient. double cimu; }; //End of TRAPPING structures Loading
