Special functions and Green’s kernels#

Special functions#

Some special functions are available from a wrapper to the Amos library provided with XLiFE++. A lot of their properties can be found on the Digital Library of Mathematical Functions (https://dlmf.nist.gov/).

Bessel and Hankel functions#

Bessel \(J_\nu\), \(Y_\nu\) and Hankel \(H^{(1)}_\nu\), \(H^{(2)}_\nu\) functions are available at any order (and also for real orders) for real or complex argument. These are solutions of differential equation:

\[\displaystyle z^2\frac{d^2w}{dz^2}+z\frac{dw}{dz}+(z^2-\nu^2)w=0.\]

and are related:

\[H^{(1)}_\nu = J\nu +i\,Y_\nu\ \mathrm{and}\ H^{(2)}_\nu = J_\nu -i\,Y_\nu.\]

  • Bessel functions of the first kind \(J_N(x)\)

    Real besselJ(Real x, Number N);
Complex besselJ(const Complex& z, Real N);
Real besselJ0(Real x);              // shortcut for order 0 and real
Complex besselJ0(const Complex& z); // shortcut for order 0 and complex
Real besselJ1(Real x);              // shortcut for order 1 and real
Complex besselJ1(const Complex& z); // shortcut for order 1 and complex
Reals besselJ0N(Real x, Number N);  // shortcut for orders 0 to N and real

  • Bessel functions of the second kind \(Y_N(x)\)

    Real besselY(Real x, Number N);
Complex besselY(const Complex& z, Real N);
Real besselY0(Real x);              // shortcut for order 0 and real
Complex besselY0(const Complex& z); // shortcut for order 0 and complex
Real besselY1(Real x);              // shortcut for order 1 and real
Complex besselY1(const Complex& z); // shortcut for order 1 and complex
Reals besselY0N(Real x, Number N);  // shortcut for orders 0 to N and real

  • Hankel functions of the first kind \(H^{(1)}_N(x)\)

    Complex hankelH1(Real x, Number N);         // general routine for real
Complex hankelH1(const Complex& z, Real N); // general routine for complex
Complex hankelH10(Real x);                  // shortcut for order 0 and real
Complex hankelH10(const Complex& z);        // shortcut for order 0 and complex e
Complex hankelH11(Real x);                  // shortcut for order 1 and real
Complex hankelH11(const Complex& z);        // shortcut for order 1 and complex
Complexes hankelH10N(Real x, Number N);     // shortcut for orders 0 to N and real

  • Hankel functions of the second kind \(H^{(2)}_N(x)\)

    Complex hankelH2(Real x, Number N);         // general routine for real
Complex hankelH2(const Complex& z, Real N); // general routine for complex
Complex hankelH20(Real x);                  // shortcut for order 0 and real
Complex hankelH20(const Complex& z);        // shortcut for order 0 and complex
Complex hankelH21(Real x);                  // shortcut for order 1 and real
Complex hankelH21(const Complex& z);        // shortcut for order 1 and complex
Complexes hankelH20N(Real x, Number N);     // shortcut for orders 0 to N and real

Modified Bessel \(I_\nu\), \(K_\nu\) functions are solutions of the differential equation:

\[\displaystyle z^2\frac{d^2w}{dz^2}+z\frac{dw}{dz}-(z^2+\nu^2)w=0.\]

  • Modified Bessel functions of the first kind \(I_N(x)\)

    Real besselI(Real x, Number N);
Complex besselI(const Complex& z, Real N);
Real besselI0(Real x);              // shortcut for order 0 and real
Complex besselI0(const Complex& z); // shortcut for order 0 and complex
Real besselI1(Real x);              // shortcut for order 1 and real
Complex besselI1(const Complex& z); // shortcut for order 1 and complex
Reals besselI0N(Real x, Number N);  // shortcut for orders 0 to N and real

  • Modified Bessel functions of the second kind \(K_N(x)\)

    Real besselK(Real x, Number N);
Complex besselK(const Complex& z, Real N);
Real besselK0(Real x);              // shortcut for order 0 and real
Complex besselK0(const Complex& z); // shortcut for order 0 and complex
Real besselK1(Real x);              // shortcut for order 1 and real
Complex besselK1(const Complex& z); // shortcut for order 1 and complex
Reals besselK0N(Real x, Number N);  // shortcut for orders 0 to N and real

  • Airy and Biry functions

    Airy \(Ai\), \(Bi\) functions are solutions of the differential equation:

    \[\displaystyle \frac{d^2w}{dz^2}=z\,w.\]
    Complex airy(Real x, DiffOpType d=_id);           // general function for real
Complex airy(const Complex& z, DiffOpType d=_id); // general function for complex
Complex biry(Real x, DiffOpType d=_id);           // general function for real
Complex biry(const Complex& z, DiffOpType d=_id); //general function for complex

Gamma and exponential integral functions#

The gamma function is defined as:

\[\displaystyle \Gamma(z)=\int_0^{+\infty}e^{-t}t^{z-1}.\]

The gamma function and related functions are available with the following routines:

Real gammaFunction(Int n);
Real gammaFunction(Real x);
Complex gammaFunction(const Complex& z);
Real logGamma(Real x);
Complex logGamma1(const Complex& z);
Complex logGamma(const Complex& z);
Real diGamma(Int n);
Real diGamma(Real x);
Complex diGamma(const Complex& z);

The exponential integrals is defined as

\[\displaystyle E_1(z)=\int_z^{+\infty}\frac{e^{-t}}{t}.\]

Related stuff is provided by the following functions:

Complex e1z(const Complex& z);  // return E1(z)
Complex eInz(const Complex& z); // return E1(z)+gamma+log(z)=\sum_{n>0} (−z)^n/nn!
Complex expzE1z(const Complex& z);  // return exp(z)∗E1(z)
Complex zexpzE1z(const Complex& z); // return z∗exp(z)∗E1(z)
Complex ascendingSeriesOfE1(const Complex& z);   // ascending series  ('small' z)
Complex continuedFractionOfE1(const Complex& z); // continued fraction ('large' z)

The erf function

\[\DeclareMathOperator{\erf}{erf} \erf z = \dfrac{2}{\sqrt {\pi}} \int_0^z e^{-t^2} dt\]

is also available from

Complex erf(Complex z);

Orthogonal polynomials#

The following orthogonal polynomials (for the inner product \(\int_{]-1,1[}w(x)\,u(x)\,v(x)\,dx\)) are available at any order

  • Legendre polynomials, orthognal with respect to the weight function \(w(x)=1\),

  • Chebyshev polynomials (\(T_n(x)=\cos(n\arccos(x))\)), orthognal with respect to the weight function \(w(x)=(1−x^2)^{-1/2}\),

  • Gegenbauer polynomials, orthognal with respect to the weight function \(w(x)=(1−x^2)^{α–1/2}\),

  • Jacobi polynomials, orthognal with respect to the weight function \(w(x)(1 − x)^\alpha(1+x)^\beta\).

void chebyshevPolynomials(Real x, Reals& val);
void gegenbauerPolynomials(Real alpha, Real x, Reals& val);
void jacobiPolynomials(Real a, Real b, Real x, Reals& val);
void jacobiPolynomials01(Real a, Real b, Real x, Reals& val); // Jacobi on [0, 1]
void legendrePolynomials(Real x, Reals& val);
void legendrePolynomialsDerivative(Real x, Reals& val);
void legendreFunctions(Real x, std::vector<Reals>& Pml);
void legendreFunctionsDerivative(Real x, const std::vector<Reals>& Plm,
                                 std::vector<Reals>& dPlm);

The number of polynomials (max order) is determined by the size of result vector (say val in exemple).

Note

Since C++17, most special functions are provided by the std library for float, double or long double types but not for complex types (see https://en.cppreference.com/w/cpp/experimental/special_functions). XLiFE++ makes no reference to these functions.

Available Green’s kernels#

XLiFE++ provides some standard Green functions as Kernel objects.

  • Laplace in free space satisfying \(-\Delta K(\mathbf{x}, \mathbf{y})=\delta_{\mathbf{x}-\mathbf{y}}\):

    \[ \begin{align}\begin{aligned}K(\mathbf{x}, \mathbf{y})=\frac{-\log(|\mathbf{x}- \mathbf{y}|)}{2\pi} \quad (2d)\\K(\mathbf{x}, \mathbf{y})=\frac{|\mathbf{x}- \mathbf{y}|}{4\pi}\quad (3d)\end{aligned}\end{align} \]

    To handle such Green’s functions:

    Kernel L2d = Laplace2dKernel();
Kernel L3d = Laplace3dKernel();

  • Helmholtz in free space satisfying \((-\Delta -k^2) K(\mathbf{x}, \mathbf{y})=\delta_{\mathbf{x}-\mathbf{y}}\):

    \[ \begin{align}\begin{aligned}\displaystyle K(\mathbf{x}, \mathbf{y})=\frac{i}{4} H^{(1)}_0(k|\mathbf{x}- \mathbf{y}|)\quad (2d)\\\displaystyle K(\mathbf{x}, \mathbf{y}) = \frac{e^{i\,k\,\left|\mathbf{x}-\mathbf{y}\right|}}{4\,\pi\,\left|\mathbf{x}-\mathbf{y}\right|} \quad (3d)\end{aligned}\end{align} \]

    where \(H^{(1)}_0\) is the Hankel function of the first kind of order 0.

    Real k=3.;
Kernel H2d = Helmholtz2dKernel(k);  // handle Hemholtz 2d kernel
Point x(0.,0.), y(1.,0.);
Complex r;
H2d(x,y,r);                         // evaluate at (x,y)

  • Helmholtz 2D in half space satisfying \((-\Delta -k^2) K(\mathbf{x}, \mathbf{y})=\delta_{\mathbf{x}-\mathbf{y}}\):

    In the 2d half space \(\mathbf{a}\mathbf{p}.\mathbf{n}>0\), where \(\mathbf{a}\) any point and \(\mathbf{n}=(-t_2,t_1)\) a normal vector to the line \(L=\{\mathbf{p}=\mathbf{a}+\alpha \mathbf{t}, \,\alpha\in\mathbb R\}\) \((\mathbf{t}=(t_1,t_2))\), the Green function

    • \(K(\mathbf{x}, \mathbf{y})=H_{2d}(\mathbf{x}, \mathbf{y}) - H_{2d}(\mathbf{x}, \mathbf{y_s})\) satisfies \(K(\mathbf{x}, \mathbf{y})=0,\, \mathbf{y}\in L\) (Dirichet condition)

    • \(K(\mathbf{x}, \mathbf{y})=H_{2d}(\mathbf{x}, \mathbf{y}) + H_{2d}(\mathbf{x}, \mathbf{y_s})\) satisfies \(\nabla_\mathbf{y}K(\mathbf{x}, \mathbf{y}).\mathbf{n}=0,\, \mathbf{y}\in L\) (Neumann condition)

    where \(H_{2d}\) is the free space Helmholtz Green function and \(\mathbf{y_s} = 2\mathbf{a} - \mathbf{yx} +2*(\mathbf{ay}.\mathbf{t})/|\mathbf{t}|^2\) the image of \(\mathbf{y}\) by symmetry. Such kernels are handled as following:

    Real k=3.;
Vector<Real> t={1.,0.};  Point a={0.,0.};
Kernel Hd = Helmholtz2dHalfPlaneKernel(k, t, a, _Dirichlet);
Kernel Hn = Helmholtz2dHalfPlaneKernel(k, t, a, _Neumann);

  • Helmhholtz 2D in strip satisfying \((-\Delta -k^2) K(\mathbf{x}, \mathbf{y})=\delta_{\mathbf{x}-\mathbf{y}}\):

    In the strip \(]-\infty,+\infty[\times\{h\}\), the Helmholtz Green function

    • \(\displaystyle K(\mathbf{x}, \mathbf{y})=\sum_{n>0}\sin\frac{n\pi x_2}{h}\sin\frac{n\pi y_2}{h}\frac{e^{-\beta_n|x_1-y_1|}}{\beta_n h}\) satisfies \(K(x,y)=0\) when \(y\) is on boundary (Dirichet condition)

    • \(\displaystyle K(\mathbf{x}, \mathbf{y})=\sum_{n\ge 0}a_n\cos\frac{n\pi x_2}{h}\cos\frac{n\pi y_2}{h}\frac{e^{-\beta_n|x_1-y_1|}}{\beta_n h}\) satisfies \(\displaystyle\frac{\partial K(\mathbf{x}, \mathbf{y})}{\partial y_2}=0\) when \(y\) is on boundary (Neumann condition)

    with \(\beta_n=\sqrt{k^2-\frac{n^2}{\pi^2h^2}}\) (\(\mathrm{Re}\beta_n\ge 0,\ \mathrm{Im}\beta_n\ge 0\)) and \(a_n=1,\ \forall n>0\), \(a_0=\frac{1}{2}\).

    When \(x_1 = y_1\), even if \(x_2 \neq y_2\), expansions have poor convergence. So, when \(x_1\) is close to \(y_1\), other expansions based on images method are used.

    Real k=3., h=2.;
Number n = 100;  // number of terms in expansion, default = 1000
Real l = 0.1;    // bound for using images method (|x1-y1| < l), default = h/10
Kernel Hs = Helmholtz2dStripKernel(_Dirichlet, k, h, n, l );

  • Harmonic regularized Maxwell 3D in free space \((\mathrm{curl}\,\mathrm{curl}-s\nabla\mathrm{div}-k^2) \mathbb{K}(\mathbf{x}, \mathbf{y})=\delta_{\mathbf{x}-\mathbf{y}}\mathbb{I}\)

    The tensor Green kernel is given by:

    \[\displaystyle \mathbb{K}(\mathbf{x}, \mathbf{y}) = H_k(\mathbf{x}, \mathbf{y})\mathbb{I}-\frac{1}{k^2}\left( \nabla^2H_k(\mathbf{x}, \mathbf{y})-\nabla^2H_{k\sqrt{s}}(\mathbf{x}, \mathbf{y})\right)\]

    with \(H_k\) the Helmholtz 3d Green’s kernel related to wave number \(k\).

    Real k=3., s=0.01;
Kernel Hs = Maxwell3dKernel(k,s);

    Warning

    This kernel is experimental. It provides only evaluation of kernel, curlx, curly, curlxy, divx, divy, divxy.

User kernels#

Users can define their own kernels, means providing the definition of all the necessary member functions: kernel, gradx, grady, gradxy, ndotgradx, ndotgrady, curlx, curly, curlxy, divx, divy, divxy, dx1, dx2, dx3.

Caution

When dealing with a matrix kernel, it may be useful to give dx1, dx2, dx3 because it is not possible to define the gradient of a matrix kernel as a Function.

To understand how it works, this is the example of Helmholtz 3D kernel:

\[K_k(\mathbf{x}, \mathbf{y}) = \frac{e^{i\,k\,\left|\mathbf{x}-\mathbf{y}\right|}}{4\,\pi\,\left|\mathbf{x}-\mathbf{y}\right|}\]
// defining Kernel K=Helmholtz3dKernel(pars) (pars containing the wave number k)
K.name="Helmholtz 3D kernel";
K.shortname="Helmz3D";
K.singularType =_r_;
K.singularOrder = -1;
K.singularCoefficient = over4pi_;
K.symmetry=_symmetric;
K.userData.push(pars);
K.kernel = Function(Helmholtz3d, K.userData);
K.gradx = Function(Helmholtz3dGradx, K.userData);
K.grady = Function(Helmholtz3dGrady, K.userData);
K.ndotgradx = Function(Helmholtz3dNxdotGradx, K.userData);
K.ndotgrady = Function(Helmholtz3dNydotGrady, K.userData);
K.singPart = new Kernel(Helmholtz3dKernelSing(K.userData));
K.regPart = new Kernel(Helmholtz3dKernelReg(K.userData));

The last 7 lines show how to define the different functions from C++ ordinary functions. Let’s see some of them in details:

  • Helmhholtz3d is the kernel function:

    // declaration
Complex Helmholtz3d(const Point& x, const Point& y, Parameters& pa = defaultParameters);
// definition (i_ and over4pi_ are available constants)
Complex Helmholtz3d(const Point& x, const Point& y, Parameters& pa)
{
  Real k = real(pa("k"));
  Real r = x.distance(y);
  Complex ikr = i_ * k * r;
  return over4pi_ * std::exp(ikr) / r;
}

  • Helmholtz3dGradx is the gradient of the kernel relative to the point x. Its analytical expression is:

    \[K_k(\mathbf{x}, \mathbf{y}) = \frac{e^{i\,k\,\left|\mathbf{x}-\mathbf{y}\right|}\,\left( i\,k\,\left|\mathbf{x}-\mathbf{y}\right| - 1\right)}{4\,\pi\,\left|\mathbf{x}-\mathbf{y}\right|^3} \left(\mathbf{x}-\mathbf{y}\right)\]
    // declaration
Vector<Complex> Helmholtz3dGradx(const Point& x, const Point& y, Parameters& pa = defaultParameters);
// definition
Vector<Complex> Helmholtz3dGradx(const Point& x, const Point& y, Parameters& pa)
{
  Complex k = pa("k");
  Real r2 = x.squareDistance(y);
  Real r = std::sqrt(r2);
  Complex ikr = i_ * k * r;
  Complex dr = (ikr - 1.) / r2;
  Vector<Complex> g1(3);
  // scaledVectorTpl computes g=s*(x-y)
  Real s = over4pi*exp(ikr) * dr / r;
  scaledVectorTpl(s, x.begin(), x.end(), y.begin(), g1.begin());
  return g1;
}

  • Let’s now talk about the kernels singPart and regPart. Both are kernel objects and have to be defined with the same mechanism from their analytical expression.

    \[K_k(\mathbf{x}, \mathbf{y}) = \frac{e^{i k \left|\mathbf{x}-\mathbf{y}\right|}}{4 \pi \left|\mathbf{x}-\mathbf{y}\right|} = \underbrace{\frac{e^{i k \left|\mathbf{x}-\mathbf{y}\right|} - 1}{4 \pi \left|\mathbf{x}-\mathbf{y}\right|}}_{\mathrm{regular\ part}} + \underbrace{\frac{1}{4 \pi \left|\mathbf{x}-\mathbf{y}\right|}}_{\mathrm{singular\ part}}\]

    Define regular part functions:

    // regular part: G_reg(k; x, y)=(exp(i*k*r)-1)/(4*pi*r)
Complex Helmholtz3dReg(const Point& x, const Point& y, Parameters& pa)
{
 Complex g;
 Real k = real(pa("k"));
 Real kr = k * x.distance(y);
 Complex ikr = Complex(0., kr);
 if (std::abs(kr) < 1.e-04)
 {
  int n=4; // for abs(kr)<1.e-4 this is a good choice for n (checked)
  g = 1 + ikr / n--;
  while (n > 1) {g = 1 + g * ikr / n--;}
  return g *= Complex(0., over4pi * k);
 }
 else return over4pi * k * (std::exp(ikr) - 1.) / kr;
}

Vector<Complex> Helmholtz3dGradxReg(const Point& x, const Point& y, Parameters& pa)
{
 Real k = real(pa("k"));
 Real r = x.distance(y);
 Complex ikr = Complex(0., k*r);
 Complex t = over4pi * (1. + std::exp(ikr)*(ikr - 1.))/ r;
 Vector<Complex> g(3);
 scaledVectorTpl(t/ r, x.begin(), x.end(), y.begin(), g.begin());
 return g;
}

Vector<Complex> Helmholtz3dGradyReg(const Point& x, const Point& y, Parameters& pa)
{
 Real k = real(pa("k"));
 Real r = x.distance(y);
 Complex ikr = Complex(0., k*r);
 Complex t = - over4pi * (1. + std::exp(ikr)*(ikr - 1.))/ r;
 Vector<Complex> g(3);
 scaledVectorTpl(t/ r, x.begin(), x.end(), y.begin(), g.begin());
 return g;
}

    Define singular part functions:

    //construct Helmholtz3d Kernel singular part: G_sing(k; x, y)=1/(4*pi*r)
Complex Helmholtz3dSing(const Point& x, const Point& y, Parameters& pa)
{
 Real r = x.distance(y);
 return over4pi/r;
}

Vector<Complex> Helmholtz3dGradxSing(const Point& x, const Point& y, Parameters& pa)
{
 Real r = x.distance(y);
 return -over4pi / (r*r);
 Complex t = -over4pi / (r*r);
 Vector<Complex> g(3);
 scaledVectorTpl(t, x.begin(), x.end(), y.begin(), g.begin());
 return g;
}

Vector<Complex> Helmholtz3dGradySing(const Point& x, const Point& y, Parameters& pa)
{
 Real r = x.distance(y);
 Complex t = over4pi / (r*r);
 Vector<Complex> g(3);
 scaledVectorTpl(t, x.begin(), x.end(), y.begin(), g.begin());
 return g;
}

Now construct Kernel objects:

Parameters pars;
pars<<Parameter(1., "k");
Kernel H3Dreg;   //regular part
H3Dreg.name="Helmholtz 3D kernel regular part";
H3Dreg.singularType =_notsingular;
H3Dreg.singularOrder = 0;
H3Dreg.singularCoefficient = over4pi;
H3Dreg.symmetry=_symmetric;
H3Dreg.userData = pars;
H3Dreg.dimPoint = 3;
H3Dreg.kernel = Function(Helmholtz3dReg, pars);
H3Dreg.gradx = Function(Helmholtz3dGradxReg, pars);
H3Dreg.grady = Function(Helmholtz3dGradyReg, pars);

Kernel H3Dsing;   //singular part
H3Dsing.name="Helmholtz 3D kernel, singular part";
H3Dsing.singularType =_r;
H3Dsing.singularOrder = -1;
H3Dsing.singularCoefficient = over4pi;
H3Dsing.symmetry=_symmetric;
H3Dsing.userData = pars;
H3Dsing.dimPoint = 3;
H3Dsing.kernel = Function(Helmholtz3dSing, pars);
H3Dsing.gradx  = Function(Helmholtz3dGradxSing, pars);
H3Dsing.grady  = Function(Helmholtz3dGradySing, pars);

Kernel H3D;   //kernel
H3D.name="Helmholtz 3D kernel";
H3D.singularType =_r;
H3D.singularOrder = -1;
H3D.singularCoefficient = over4pi;
H3D.symmetry=_symmetric;
H3D.userData = pars;
H3D.dimPoint = 3;
H3D.kernel = Function(Helmholtz3d, pars);
H3D.gradx  = Function(Helmholtz3dGradx, pars);
H3D.grady  = Function(Helmholtz3dGrady, pars);
H3D.regPart  = &H3Dreg;
H3D.singPart = &H3Dsing;

Warning

If singular and regular part of kernels are not defined, some computations will not be available.

Tip

New kernels developped by users may be of interest for others. So contact administrators, they will be happy to integrate new works in XLiFE++.

PDE exact solutions#

The field scattered by a disk of radius \(R\) of an incoming field \(\Phi_w=e^{ikx}\) with Neumann or Dirichlet boundary conditions on the boundary are given by the Mie series expansions:

\[\begin{split}\displaystyle u(r,\theta)=-\sum_{n\in\mathbb{Z}}(-i)^n\frac{J_n(kR)}{H^{(1)}_n(kR)}H^{(1)}_n(kr)e^{in\theta}\quad \mathrm{(Dirichlet)}\\ \displaystyle u(r,\theta)=-\sum_{n\in\mathbb{Z}}(-i)^n\frac{J_n'(kR)}{H^{(1)'}_n(kR)}H^{(1)}_n(kr)e^{in\theta}\quad \mathrm{(Neumann)}\\\end{split}\]

These solutions are available in XLiFE++:

Real k=10., R=1.;  // wave number and disk radius
Number nmax=30;    // number of terms in expansion (2*nmax+1)
Parameters pars;
pars<<Parameter(k,"k")<<Parameter(R,"radius")<<Parameter(nmax,"nmax");
Point P(1.,1.);
Complex scatteredFieldDiskDirichlet(P,pars);
Complex scatteredFieldDiskNeumann(P,pars);

If parameters are not set, the default values are k=1, R=1 and nmax=50.

The field scattered by a solid sphere of radius \(R\) of an incoming field \(\Phi_w=e^{ikx}\) with Neumann or Dirichlet boundary conditions on the sphere may be obtained as series expansions involving spherical Bessel functions and Legendre functions (spherical harmonics). XLiFE++ provides the two functions:

Complex scatteredFieldSphereNeumann(const Point& x, Parameters& param);
Complex scatteredFieldSphereDirichlet(const Point& x, Parameters& param);

The parameters k,r,nmax are the same as the previous ones.

The field scattered by a solid sphere of radius R of an incoming field \(\Phi_w=e^{ikz}\) with tangential boundary condition \(E\times n=0\) (electric field) on the sphere is also available with the following function:

Complexes scatteredFieldMaxwellExn(const Point& x, Parameters& param);

At last, the spherical harmonics \(Y_l^m\) up to order \(n\) on the unit sphere and their gradient are available with the following functions:

void sphericalHarmonics(const Point& x, vector<Complexes>& Yml);
void sphericalHarmonicsSurfaceGrad(const Point& x, vector<vector<Complexes> >& gradylm);