GeographicLib  1.48
SphericalEngine.cpp File Reference

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

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

Go to the source code of this file.

## Namespaces

GeographicLib
Namespace for GeographicLib.

## Detailed Description

Implementation for GeographicLib::SphericalEngine class.

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.