Implementation for GeographicLib::SphericalEngine class. More...

#include <GeographicLib/SphericalEngine.hpp>
#include <GeographicLib/CircularEngine.hpp>
#include <GeographicLib/Utility.hpp>

Namespaces
namespace	GeographicLib
	Namespace for GeographicLib.

Detailed Description

Implementation for GeographicLib::SphericalEngine class.

Copyright (c) Charles Karney (2011-2022) charl.nosp@m.es@k.nosp@m.arney.nosp@m..com and licensed under the MIT/X11 License. For more information, see https://geographiclib.sourceforge.io/

The general sum is

 V(r, theta, lambda) = sum(n = 0..N) sum(m = 0..n)
   q^(n+1) * (C[n,m] * cos(m*lambda) + S[n,m] * sin(m*lambda)) * P[n,m](t)

where t = cos(theta), q = a/r. In addition write u = sin(theta).

P[n,m] is a normalized associated Legendre function of degree n and order m. Here the formulas are given for full normalized functions (usually denoted Pbar).

Rewrite outer sum

 V(r, theta, lambda) = sum(m = 0..N) * P[m,m](t) * q^(m+1) *
    [Sc[m] * cos(m*lambda) + Ss[m] * sin(m*lambda)]

where the inner sums are

   Sc[m] = sum(n = m..N) q^(n-m) * C[n,m] * P[n,m](t)/P[m,m](t)
   Ss[m] = sum(n = m..N) q^(n-m) * S[n,m] * P[n,m](t)/P[m,m](t)

Evaluate sums via Clenshaw method. The overall framework is similar to Deakin with the following changes:

Clenshaw summation is used to roll the computation of cos(m*lambda) and sin(m*lambda) into the evaluation of the outer sum (rather than independently computing an array of these trigonometric terms).
Scale the coefficients to guard against overflow when N is large.

For the general framework of Clenshaw, see http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html

Let

    S = sum(k = 0..N) c[k] * F[k](x)
    F[n+1](x) = alpha[n](x) * F[n](x) + beta[n](x) * F[n-1](x)

Evaluate S with

    y[N+2] = y[N+1] = 0
    y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]
    S = c[0] * F[0] + y[1] * F[1] + beta[1] * F[0] * y[2]

IF F[0](x) = 1 and beta(0,x) = 0, then F[1](x) = alpha(0,x) and we can continue the recursion for y[k] until y[0], giving

    S = y[0]

Evaluating the inner sum

 l = n-m; n = l+m
 Sc[m] = sum(l = 0..N-m) C[l+m,m] * q^l * P[l+m,m](t)/P[m,m](t)
 F[l] = q^l * P[l+m,m](t)/P[m,m](t)

Holmes + Featherstone, Eq. (11), give

   P[n,m] = sqrt((2*n-1)*(2*n+1)/((n-m)*(n+m))) * t * P[n-1,m] -
            sqrt((2*n+1)*(n+m-1)*(n-m-1)/((n-m)*(n+m)*(2*n-3))) * P[n-2,m]

thus

   alpha[l] = t * q * sqrt(((2*n+1)*(2*n+3))/
                           ((n-m+1)*(n+m+1)))
   beta[l+1] = - q^2 * sqrt(((n-m+1)*(n+m+1)*(2*n+5))/
                            ((n-m+2)*(n+m+2)*(2*n+1)))

In this case, F[0] = 1 and beta[0] = 0, so the Sc[m] = y[0].

Evaluating the outer sum

 V = sum(m = 0..N) Sc[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
   + sum(m = 0..N) Ss[m] * q^(m+1) * cos(m*lambda) * P[m,m](t)
 F[m] = q^(m+1) * cos(m*lambda) * P[m,m](t) [or sin(m*lambda)]

Holmes + Featherstone, Eq. (13), give

   P[m,m] = u * sqrt((2*m+1)/((m>1?2:1)*m)) * P[m-1,m-1]

also, we have

   cos((m+1)*lambda) = 2*cos(lambda)*cos(m*lambda) - cos((m-1)*lambda)

thus

   alpha[m] = 2*cos(lambda) * sqrt((2*m+3)/(2*(m+1))) * u * q
            =   cos(lambda) * sqrt( 2*(2*m+3)/(m+1) ) * u * q
   beta[m+1] = -sqrt((2*m+3)*(2*m+5)/(4*(m+1)*(m+2))) * u^2 * q^2
               * (m == 0 ? sqrt(2) : 1)

Thus

 F[0] = q                                [or 0]
 F[1] = cos(lambda) * sqrt(3) * u * q^2  [or sin(lambda)]
 beta[1] = - sqrt(15/4) * u^2 * q^2

Here is how the various components of the gradient are computed

Differentiate wrt r

   d q^(n+1) / dr = (-1/r) * (n+1) * q^(n+1)

so multiply C[n,m] by n+1 in inner sum and multiply the sum by -1/r.

Differentiate wrt lambda

   d cos(m*lambda) = -m * sin(m*lambda)
   d sin(m*lambda) =  m * cos(m*lambda)

so multiply terms by m in outer sum and swap sine and cosine variables.

Differentiate wrt theta

  dV/dtheta = V' = -u * dV/dt = -u * V'

here ' denotes differentiation wrt theta.

   d/dtheta (Sc[m] * P[m,m](t)) = Sc'[m] * P[m,m](t) + Sc[m] * P'[m,m](t)

Now P[m,m](t) = const * u^m, so P'[m,m](t) = m * t/u * P[m,m](t), thus

   d/dtheta (Sc[m] * P[m,m](t)) = (Sc'[m] + m * t/u * Sc[m]) * P[m,m](t)

Clenshaw recursion for Sc[m] reads

    y[k] = alpha[k] * y[k+1] + beta[k+1] * y[k+2] + c[k]

Substituting alpha[k] = const * t, alpha'[k] = -u/t * alpha[k], beta'[k] = c'[k] = 0 gives

    y'[k] = alpha[k] * y'[k+1] + beta[k+1] * y'[k+2] - u/t * alpha[k] * y[k+1]

Finally, given the derivatives of V, we can compute the components of the gradient in spherical coordinates and transform the result into cartesian coordinates.

Definition in file SphericalEngine.cpp.

Namespaces

Detailed Description